PIRE analyses for paper.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as sf\n",
    "import plotly.graph_objects as go\n",
    "from IPython.display import Image\n",
    "from pi_for_temperature_notebook import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Veg6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data\n",
    "veg6_data = pd.read_csv(\"variant-calls_vegprop-expt_2020-04-07_tables.csv\",header=0,index_col=0)\n",
    "veg6_options = get_veg_prop_args()\n",
    "\n",
    "#filter SNPs\n",
    "filtered_veg6_data = pd.concat([veg6_data.loc[veg6_data.freqPropReplicateA >= 0.03],\n",
    "                           veg6_data.loc[veg6_data.freqPropReplicateB >= 0.03]]).drop_duplicates()\n",
    "\n",
    "filtered_veg6_data['DP'] = filtered_veg6_data['ADReplicateA'] + filtered_veg6_data['ADReplicateB'] + filtered_veg6_data['RDReplicateA'] + filtered_veg6_data['RDReplicateB']\n",
    "filtered_veg6_data['AD'] = filtered_veg6_data['ADReplicateA'] + filtered_veg6_data['ADReplicateB']\n",
    "\n",
    "#calculate frequency of aggregated read counts\n",
    "filtered_veg6_data['freqProp'] = filtered_veg6_data['AD'] / filtered_veg6_data['DP']\n",
    "\n",
    "\n",
    "#lineage: set to 0 then assign each unique lineage to an unique integer\n",
    "filtered_veg6_data['lineage_factor'] = 0\n",
    "j = 0\n",
    "for lineage in filtered_veg6_data.lineage.unique():\n",
    "    filtered_veg6_data.loc[filtered_veg6_data.lineage == lineage,['lineage_factor']] = j\n",
    "    j += 1\n",
    "\n",
    "#species: Set ACMV = 0 and EACMCV = 1\n",
    "filtered_veg6_data['species'] = 0\n",
    "filtered_veg6_data.loc[filtered_veg6_data.chrom == 'EACMCV DNA-A',['species']] = 1\n",
    "filtered_veg6_data.loc[filtered_veg6_data.chrom == 'EACMCV DNA-B',['species']] = 1\n",
    "\n",
    "#segment: Set DNA-A = 0 and DNA-B = 1\n",
    "filtered_veg6_data['segment'] = 0\n",
    "filtered_veg6_data.loc[filtered_veg6_data.chrom == 'ACMV DNA-B',['segment']] = 1\n",
    "filtered_veg6_data.loc[filtered_veg6_data.chrom == 'EACMCV DNA-B',['segment']] = 1\n",
    "\n",
    "#SEGs treatment: Virus Only:0, SEGS-1:1,SEGS-2:2\n",
    "filtered_veg6_data['SEGs_treatment'] = 0\n",
    "filtered_veg6_data.loc[filtered_veg6_data.segsTreatment2 == 'SEGS-1','SEGs_treatment'] = 1\n",
    "filtered_veg6_data.loc[filtered_veg6_data.segsTreatment2 == 'SEGS-2','SEGs_treatment'] = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg6_options['frequency'] = 'freqPropMeanNoNA'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Erik\\OneDrive\\DuffyLab\\PIRE\\pi_for_temperature_notebook.py:158: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  site_data[frequency] = [x/dp for x in site_data['AD']]\n"
     ]
    }
   ],
   "source": [
    "veg6_avg = average_groups(filtered_veg6_data,['species','segment','pos','alt','passage','ref','lineage_factor','SEGs_treatment'])\n",
    "\n",
    "veg6_avg.reset_index(inplace=True)\n",
    "veg6_avg['pos'] = veg6_avg.pos.astype(int)\n",
    "veg6_pi_df = get_group_pis(veg6_avg,\n",
    "                     options=veg6_options,\n",
    "                     group_by=['species','segment','passage','lineage_factor','SEGs_treatment'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Least squares summary:\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                     pi   R-squared:                       0.179\n",
      "Model:                            OLS   Adj. R-squared:                  0.139\n",
      "Method:                 Least Squares   F-statistic:                     3.077\n",
      "Date:                Mon, 10 Aug 2020   Prob (F-statistic):           0.000114\n",
      "Time:                        12:14:16   Log-Likelihood:                -851.96\n",
      "No. Observations:                 323   AIC:                             1736.\n",
      "Df Residuals:                     307   BIC:                             1796.\n",
      "Df Model:                          15                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "=================================================================================================\n",
      "                                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------------\n",
      "Intercept                        -0.0055      0.735     -0.007      0.994      -1.447       1.436\n",
      "C(SEGs_treatment)[T.1]            0.3150      0.354      0.889      0.374      -0.379       1.009\n",
      "C(SEGs_treatment)[T.2]            1.0696      0.494      2.164      0.030       0.101       2.038\n",
      "C(species)[T.1]                   0.3634      0.574      0.633      0.527      -0.762       1.489\n",
      "C(lineage_factor)[T.1]            1.4820      0.726      2.041      0.041       0.059       2.905\n",
      "C(lineage_factor)[T.2]            0.6129      0.684      0.896      0.370      -0.727       1.953\n",
      "C(lineage_factor)[T.3]            0.6416      0.663      0.968      0.333      -0.657       1.940\n",
      "C(lineage_factor)[T.4]           -0.6529      0.325     -2.010      0.044      -1.290      -0.016\n",
      "C(lineage_factor)[T.5]           -0.1497      0.251     -0.596      0.551      -0.642       0.342\n",
      "C(lineage_factor)[T.6]            0.3912      0.340      1.150      0.250      -0.275       1.058\n",
      "C(lineage_factor)[T.7]            0.7264      0.473      1.535      0.125      -0.201       1.654\n",
      "C(lineage_factor)[T.8]           -0.9060      0.497     -1.824      0.068      -1.879       0.067\n",
      "C(lineage_factor)[T.9]            1.8955      1.246      1.521      0.128      -0.547       4.338\n",
      "C(lineage_factor)[T.10]          -0.8104      0.558     -1.452      0.147      -1.905       0.284\n",
      "C(lineage_factor)[T.11]           0.8905      0.805      1.107      0.268      -0.687       2.468\n",
      "C(species)[0]:C(segment)[T.1]     0.0131      0.453      0.029      0.977      -0.875       0.901\n",
      "C(species)[1]:C(segment)[T.1]    -0.3888      0.617     -0.630      0.529      -1.598       0.820\n",
      "passage                           0.6564      0.120      5.459      0.000       0.421       0.892\n",
      "==============================================================================\n",
      "Omnibus:                      276.848   Durbin-Watson:                   2.003\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             7002.700\n",
      "Skew:                           3.396   Prob(JB):                         0.00\n",
      "Kurtosis:                      24.776   Cond. No.                     7.15e+15\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "[2] The smallest eigenvalue is 1.33e-28. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n",
      "\n",
      "Anova table:\n",
      "                          df       sum_sq     mean_sq          F        PR(>F)\n",
      "C(SEGs_treatment)        2.0    49.292435   24.646218   2.047068  1.308716e-01\n",
      "C(species)               1.0     3.097202    3.097202   0.257248  6.123815e-01\n",
      "C(lineage_factor)       11.0   211.682384   19.243853   1.598358  9.806364e-02\n",
      "C(species):C(segment)    2.0     7.994833    3.997417   0.332018  7.177318e-01\n",
      "passage                  1.0   538.145622  538.145622  44.697351  1.085360e-10\n",
      "Residual               307.0  3696.208032   12.039766        NaN           NaN\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\base\\model.py:1832: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 17, but rank is 15\n",
      "  'rank is %d' % (J, J_), ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "lm = sf.ols('pi ~ passage + C(SEGs_treatment) + (C(species)/C(segment)) + C(lineage_factor)',data=veg6_pi_df).fit(cov_type='HC1')\n",
    "print(\"Least squares summary:\")\n",
    "print(lm.summary())\n",
    "\n",
    "print(\"\\nAnova table:\")\n",
    "veg6_table = sm.stats.anova_lm(lm)\n",
    "print(veg6_table)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg6_table.to_csv(\"Veg6_full.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Erik\\OneDrive\\DuffyLab\\PIRE\\pi_for_temperature_notebook.py:158: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  site_data[frequency] = [x/dp for x in site_data['AD']]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Least squares summary:\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                     pi   R-squared:                       0.576\n",
      "Model:                            OLS   Adj. R-squared:                  0.334\n",
      "Method:                 Least Squares   F-statistic:                     3.162\n",
      "Date:                Mon, 10 Aug 2020   Prob (F-statistic):             0.0878\n",
      "Time:                        12:36:06   Log-Likelihood:                -42.683\n",
      "No. Observations:                  12   AIC:                             95.37\n",
      "Df Residuals:                       7   BIC:                             97.79\n",
      "Df Model:                           4                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "=================================================================================================\n",
      "                                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------------\n",
      "Intercept                        15.5757     11.317      1.376      0.169      -6.605      37.756\n",
      "C(species)[T.1]                  -2.2461      9.425     -0.238      0.812     -20.719      16.226\n",
      "C(species)[0]:C(segment)[T.1]   -17.0654     13.116     -1.301      0.193     -42.773       8.642\n",
      "C(species)[1]:C(segment)[T.1]     5.3539      5.379      0.995      0.320      -5.188      15.896\n",
      "passage                           3.3823      1.405      2.408      0.016       0.629       6.135\n",
      "==============================================================================\n",
      "Omnibus:                        0.560   Durbin-Watson:                   2.818\n",
      "Prob(Omnibus):                  0.756   Jarque-Bera (JB):                0.590\n",
      "Skew:                           0.315   Prob(JB):                        0.744\n",
      "Kurtosis:                       2.115   Cond. No.                         18.0\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "\n",
      "Anova table:\n",
      "                        df      sum_sq     mean_sq         F    PR(>F)\n",
      "C(species)             1.0   52.857880   52.857880  0.428516  0.533628\n",
      "C(species):C(segment)  2.0  441.499625  220.749813  1.789606  0.235644\n",
      "passage                1.0  679.724714  679.724714  5.510489  0.051280\n",
      "Residual               7.0  863.457465  123.351066       NaN       NaN\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\stats.py:1450: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=12\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "veg6_avg = average_groups(filtered_veg6_data,['species','segment','pos','alt','passage','ref'])\n",
    "\n",
    "veg6_avg.reset_index(inplace=True)\n",
    "veg6_avg['pos'] = veg6_avg.pos.astype(int)\n",
    "veg6_pi_df = get_group_pis(veg6_avg,\n",
    "                     options=veg6_options,\n",
    "                     group_by=['species','segment','passage'])\n",
    "\n",
    "lm = sf.ols('pi ~ passage + C(species)/C(segment)',data=veg6_pi_df).fit(cov_type='HC1')\n",
    "print(\"Least squares summary:\")\n",
    "print(lm.summary())\n",
    "\n",
    "print(\"\\nAnova table:\")\n",
    "veg6_table = sm.stats.anova_lm(lm)\n",
    "print(veg6_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg6_table.to_csv(\"Veg6_minimal.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Veg2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py:3063: DtypeWarning: Columns (9,58) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n"
     ]
    }
   ],
   "source": [
    "data = pd.read_csv(\"variant-calls_temperature-expt_2020-04-07_tables.csv\", header=0, index_col=0)\n",
    "veg2_options = get_veg_prop_args()\n",
    "\n",
    "data_no_32 = data.loc[data.temperature != 32]\n",
    "data_no_32_both_present = data_no_32.dropna(subset=['ADReplicateA','ADReplicateB'])\n",
    "filtered_data = pd.concat([data_no_32_both_present.loc[data_no_32_both_present.freqPropReplicateA >= 0.3],\n",
    "                           data_no_32_both_present.loc[data_no_32_both_present.freqPropReplicateB >= 0.03]]).drop_duplicates()\n",
    "\n",
    "filtered_data['DP'] = filtered_data['ADReplicateA'] + filtered_data['ADReplicateB'] + filtered_data['RDReplicateA'] + filtered_data['RDReplicateB']\n",
    "filtered_data['AD'] = filtered_data['ADReplicateA'] + filtered_data['ADReplicateB']\n",
    "\n",
    "filtered_data['lineage_factor'] = 0\n",
    "\n",
    "j = 0\n",
    "for lineage in filtered_data.lineage.unique():\n",
    "    filtered_data.loc[filtered_data.lineage == lineage,['lineage_factor']] = j\n",
    "    j += 1\n",
    "\n",
    "#species: Set ACMV = 0 and EACMCV = 1\n",
    "filtered_data['species'] = 0\n",
    "filtered_data.loc[filtered_data.chrom == 'EACMCV DNA-A',['species']] = 1\n",
    "filtered_data.loc[filtered_data.chrom == 'EACMCV DNA-B',['species']] = 1\n",
    "\n",
    "#segment: Set DNA-A = 0 and DNA-B = 1\n",
    "filtered_data['segment'] = 0\n",
    "filtered_data.loc[filtered_data.chrom == 'ACMV DNA-B',['segment']] = 1\n",
    "filtered_data.loc[filtered_data.chrom == 'EACMCV DNA-B',['segment']] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg2_options['frequency'] = 'freqPropMeanNoNA'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg2_avg = average_groups(filtered_data,['species','segment','pos','alt','passage','ref','lineage_factor','temperature'])\n",
    "\n",
    "veg2_avg.reset_index(inplace=True)\n",
    "veg2_avg['pos'] = veg2_avg.pos.astype(int)\n",
    "\n",
    "veg2_pi_df = get_group_pis(veg2_avg,\n",
    "                     options=veg2_options,\n",
    "                     group_by=['species','segment','passage','lineage_factor','temperature'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Least squares summary:\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                     pi   R-squared:                       0.400\n",
      "Model:                            OLS   Adj. R-squared:                  0.306\n",
      "Method:                 Least Squares   F-statistic:                     5.416\n",
      "Date:                Mon, 10 Aug 2020   Prob (F-statistic):           2.38e-05\n",
      "Time:                        12:15:40   Log-Likelihood:                -123.36\n",
      "No. Observations:                  67   AIC:                             266.7\n",
      "Df Residuals:                      57   BIC:                             288.8\n",
      "Df Model:                           9                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "=================================================================================================\n",
      "                                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------------\n",
      "Intercept                         0.2489      0.664      0.375      0.708      -1.053       1.551\n",
      "C(temperature)[T.30]              0.0719      0.347      0.207      0.836      -0.608       0.752\n",
      "C(species)[T.1]                  -0.6877      0.691     -0.996      0.319      -2.041       0.666\n",
      "C(lineage_factor)[T.1]            0.1715      0.450      0.381      0.703      -0.711       1.054\n",
      "C(lineage_factor)[T.2]           -0.3143      0.439     -0.716      0.474      -1.175       0.546\n",
      "C(lineage_factor)[T.3]            0.7574      0.747      1.014      0.311      -0.707       2.222\n",
      "C(lineage_factor)[T.4]           -0.2760      0.504     -0.547      0.584      -1.264       0.712\n",
      "C(lineage_factor)[T.5]           -0.4095      0.359     -1.140      0.254      -1.114       0.295\n",
      "C(species)[0]:C(segment)[T.1]    -0.9141      0.599     -1.527      0.127      -2.088       0.259\n",
      "C(species)[1]:C(segment)[T.1]     0.2448      0.559      0.438      0.662      -0.852       1.341\n",
      "passage                           1.4147      0.271      5.220      0.000       0.883       1.946\n",
      "==============================================================================\n",
      "Omnibus:                       29.208   Durbin-Watson:                   2.187\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               67.413\n",
      "Skew:                           1.383   Prob(JB):                     2.30e-15\n",
      "Kurtosis:                       7.061   Cond. No.                     3.44e+16\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "[2] The smallest eigenvalue is 3.7e-31. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n",
      "\n",
      "Anova table:\n",
      "                         df      sum_sq    mean_sq          F    PR(>F)\n",
      "C(temperature)          1.0    2.283153   2.283153   0.834818  0.364732\n",
      "C(species)              1.0    0.248183   0.248183   0.090746  0.764327\n",
      "C(lineage_factor)       5.0   17.651142   3.530228   1.290802  0.280653\n",
      "C(species):C(segment)   2.0   10.101245   5.050622   1.846722  0.167068\n",
      "passage                 1.0   76.621713  76.621713  28.016158  0.000002\n",
      "Residual               57.0  155.889954   2.734911        NaN       NaN\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\base\\model.py:1832: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 10, but rank is 9\n",
      "  'rank is %d' % (J, J_), ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "lm = sf.ols('pi ~ passage + C(temperature) + (C(species)/C(segment)) + C(lineage_factor)',data=veg2_pi_df).fit(cov_type='HC1')\n",
    "print(\"Least squares summary:\")\n",
    "print(lm.summary())\n",
    "\n",
    "print(\"\\nAnova table:\")\n",
    "veg2_table = sm.stats.anova_lm(lm)\n",
    "print(veg2_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg2_table.to_csv(\"Veg2_full.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Erik\\OneDrive\\DuffyLab\\PIRE\\pi_for_temperature_notebook.py:158: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  site_data[frequency] = [x/dp for x in site_data['AD']]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Least squares summary:\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                     pi   R-squared:                       0.834\n",
      "Model:                            OLS   Adj. R-squared:                  0.740\n",
      "Method:                 Least Squares   F-statistic:                     16.42\n",
      "Date:                Mon, 10 Aug 2020   Prob (F-statistic):            0.00115\n",
      "Time:                        12:15:55   Log-Likelihood:                -31.674\n",
      "No. Observations:                  12   AIC:                             73.35\n",
      "Df Residuals:                       7   BIC:                             75.77\n",
      "Df Model:                           4                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "=================================================================================================\n",
      "                                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------------\n",
      "Intercept                        -1.4767      3.978     -0.371      0.711      -9.274       6.321\n",
      "C(species)[T.1]                  -4.3623      2.749     -1.587      0.113      -9.751       1.026\n",
      "C(species)[0]:C(segment)[T.1]    -4.0079      4.432     -0.904      0.366     -12.695       4.679\n",
      "C(species)[1]:C(segment)[T.1]     1.4958      2.571      0.582      0.561      -3.544       6.535\n",
      "passage                           9.0743      1.393      6.513      0.000       6.344      11.805\n",
      "==============================================================================\n",
      "Omnibus:                        0.115   Durbin-Watson:                   2.884\n",
      "Prob(Omnibus):                  0.944   Jarque-Bera (JB):                0.320\n",
      "Skew:                           0.130   Prob(JB):                        0.852\n",
      "Kurtosis:                       2.244   Cond. No.                         10.2\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "\n",
      "Anova table:\n",
      "                        df      sum_sq     mean_sq          F    PR(>F)\n",
      "C(species)             1.0    7.780825    7.780825   0.395130  0.549568\n",
      "C(species):C(segment)  2.0   27.450891   13.725446   0.697013  0.529599\n",
      "passage                1.0  658.749060  658.749060  33.452969  0.000675\n",
      "Residual               7.0  137.842577   19.691797        NaN       NaN\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\stats.py:1450: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=12\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "veg2_avg = average_groups(filtered_data,['species','segment','pos','alt','passage','ref'])\n",
    "\n",
    "veg2_avg.reset_index(inplace=True)\n",
    "veg2_avg['pos'] = veg2_avg.pos.astype(int)\n",
    "\n",
    "veg2_pi_df = get_group_pis(veg2_avg,\n",
    "                     options=veg2_options,\n",
    "                     group_by=['species','segment','passage'])\n",
    "\n",
    "lm = sf.ols('pi ~ passage + (C(species)/C(segment))',data=veg2_pi_df).fit(cov_type='HC1')\n",
    "print(\"Least squares summary:\")\n",
    "print(lm.summary())\n",
    "\n",
    "print(\"\\nAnova table:\")\n",
    "veg2_table = sm.stats.anova_lm(lm)\n",
    "print(veg2_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg2_table.to_csv(\"Veg2_minimal.csv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}