import statsmodels.api as sm
from statsmodels.formula.api import ols
import scipy.stats as st
import pandas as pd
import numpy as np

dfl = pd.DataFrame({'feed':['A', 'B','C','D']*4,
          'crop': [60, 62, 68, 56, 59, 67, 67, 60,
                    59, 63, 71, 59, 58, 64, 70, 57]})
dfw = pd.DataFrame({'A':[60,59,59,58], 'B':[62,67,63,64],
                    'C':[68,67,71,70], 'D':[56,60,59,57]})

#---------------- scipy.stats --------------------------
F1, p1 = st.f_oneway(dfw['A'], dfw['B'], dfw['C'], dfw['D'])
print(F1, p1)
#---------------- statsmodels ols ----------------------
lm1 = ols('crop~C(feed)', data=dfl)
res1=lm1.fit()
aov_table1 = sm.stats.anova_lm(res1, type=3)
print(aov_table1)
print(res1.summary())

dfl.boxplot('crop', by='feed', figsize=(12, 8))

#############################################################

x = np.array([0.005, 0.600, 1.325, 1.505, 1.520, 2.000, 2.255, 2.575, 
              3.000, 3.300, 3.500] )
y = np.array([0.005, 0.005, 0.255, 0.505, 0.525, 1.000, 1.250, 1.505, 
              1.750, 2.750])

F2, p2 = st.f_oneway(x, y)
print(F2, p2)

group = np.append(np.zeros(len(x)), np.ones(len(y))) 
df2 = pd.DataFrame({'x': group, 'dep': np.append(x,y)})
res2 = ols('dep~C(group)', data=df2).fit()
aov_table2 = sm.stats.anova_lm(res2, type=3)
print(aov_table2)