# -*- coding: utf-8 -*-
"""
Created on Wed Dec  1 16:57:31 2021

@author: Sim
"""
# import naginterfaces.library.anova as aov
import scipy.stats as st
import pandas as pd
import numpy as np

def contrast(dfl, ctrt, indep):  # assumes long format data 
  mm = dfl.groupby([indep]).mean()  # group mean
  vv = dfl.groupby([indep]).var(ddof=1) # group std
  nn = dfl.groupby([indep]).count()  # group nobs
  N = np.sum(nn)       # total nobs
  dfE = N-len(mm)      # 자유도
  SSE = np.sum((nn-1)*vv) 
  MSE = SSE/dfE
  statist = np.sum(mm*ctrt) / np.sqrt(MSE*np.sum(ctrt**2/nn))
  pvalue = (1-st.t.cdf(np.abs(statist), dfE))*2 # for two-sided test
  return([statist[0], pvalue[0]])

#-----------------  
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]})
ctrt = np.array([3,-1,-1,-1]).reshape(-1,1)  
# groupby에 의해 생성된 값이 열벡터로 내부계산에서 차원 일치를 위해 열로
t1, p1 = contrast(dfl, ctrt, 'feed')