# -*- coding: utf-8 -*-
"""
Created on Thu Dec  2 11:14:42 2021

@author: Sim
"""
import scipy.stats as st
import pandas as pd
import numpy as np

def contrast2(dfl, ctrt, indep, alternative='two-sided', alpha=0.05):  
# dfl: long format 그룹변수 indep와 반응변수만 포함한 nx2 차원 DataFrame. 
  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
  
  se = np.sqrt(MSE*np.sum(ctrt**2/nn))
  cntrt = np.sum(mm*ctrt)
  lci = cntrt - st.t.ppf(1-alpha/2, dfE)*se
  uci = cntrt + st.t.ppf(1-alpha/2, dfE)*se
  
  statist = cntrt / se
  if alternative == 'greater':
    pvalue = 1-st.t.cdf(np.abs(statist), dfE)
  elif alternative == 'less':
    pvalue = st.t.cdf(np.abs(statist), dfE)
  else:
    pvalue = (1-st.t.cdf(np.abs(statist), dfE))*2 
  return([statist[0], pvalue[0], (lci[0], uci[0])])
#  return([statist, pvalue, (lci, uci)])
#-----------------  

if __name__ == '__main__':
  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에 의해 생성된 값이 열벡터로 내부계산에서 차원 일치를 위해 열로
  t, p, ci1 = contrast2(dfl, ctrt, 'feed')
  print(ci1)
  t, p, ci2 = contrast2(dfl, ctrt, 'feed', alpha=0.01)
  print(ci2)