8.10.1 Bartletts test i gruppespecifik regression
xxxxxxxxxx
import numpy as np
import pandas as pd
import scipy.stats as st
import matplotlib.pyplot as plt
from statsmodels.formula.api import ols
# funktionen bartlettList defineres
def bartlettList(lmUDlist:"liste_olsOutputs"):
# Bartletts test baseret på liste med ols outputs
k=len(lmUDlist)
if k>2:
s2=np.zeros(k); df=np.zeros(k)
for j in range(0,k):
s2[j]=lmUDlist[j].mse_resid
df[j]=lmUDlist[j].df_resid
df0=sum(df); s20=sum(df*s2)/df0
Ba=(df0*np.log(s20)-sum(df*np.log(s2)))/(1+(sum(1/df)-1/df0)/(3*(k-1)))
pval=1-st.chi2.cdf(Ba,k-1)
return(pd.DataFrame(np.array([Ba,k-1,pval]).reshape(1,-1),
columns=['Statistic','DF','Pvalue'],index=['']))
else:
s21=lmUDlist[0].mse_resid; s22=lmUDlist[1].mse_resid
df1=lmUDlist[0].df_resid; df2=lmUDlist[1].df_resid
fstat=s21/s22
pval=2*np.min([st.f.cdf(fstat,df1,df2),1-st.f.cdf(fstat,df1,df2)])
return(pd.DataFrame(np.array([fstat,df1,df2,pval]).reshape(1,-1),
columns=['Statistic','DF1','DF2','Pvalue'],index=['']))
# data indskrives
haemeff0=np.array([75.68,79.58,83.89,85.55,88.11,76.49,84.76,85.58,
86.26,91.43,77.38,86.92,88.85,89.92,93.00])
logKonc0=np.log(np.tile(np.array([0.1,0.25,0.50,0.75,1.00]),3))
gr0=np.repeat(["INH1","INH2","INH3"],5)
# datatabel dannes
# (hvis data indlæses med read_csv har man allerede en datatabel)
datahaemmer=pd.DataFrame({'haemeff':haemeff0,'logKonc':logKonc0,'gr':gr0})
print("Opstart er gennemført: haemeff, logKonc og gr er indskrevet i tabel datahaemmer")
xxxxxxxxxx
# Opstart ovenfor skal være kørt. Data er i datatabellen datahaemmer
# Figur med data
plt.plot(logKonc0[gr0=='INH1'],haemeff0[gr0=='INH1'],'o')
plt.plot(logKonc0[gr0=='INH2'],haemeff0[gr0=='INH2'],'ro')
plt.plot(logKonc0[gr0=='INH3'],haemeff0[gr0=='INH3'],'go')
# gruppespecifik model analyseres
lmUD1=ols(data=datahaemmer[datahaemmer.G=="INH1"],formula='haemeff~logKonc').fit()
lmUD2=ols(data=datahaemmer[datahaemmer.G=="INH2"],formula='haemeff~logKonc').fit()
lmUD3=ols(data=datahaemmer[datahaemmer.G=="INH3"],formula='haemeff~logKonc').fit()
# bartletts test laves
print(bartlettList([lmUD1,lmUD2,lmUD3]))
# regressionslinjer indtegnes
par=lmUD1.params
plt.axline([0.5,par[0]+0.5*par[1]],slope=par[1])
par=lmUD2.params
plt.axline([0.5,par[0]+0.5*par[1]],slope=par[1],color='r')
par=lmUD3.params
plt.axline([0.5,par[0]+0.5*par[1]],slope=par[1],color='g')
plt.show()
8.10.2 F-tests og parametertabel
xxxxxxxxxx
import numpy as np
import pandas as pd
import scipy.stats as st
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
# funktionen summaryLM defineres
def summaryLM(lmUD:"ols_output"):
# lave summary af estimeret model i ols-output
pd.options.display.float_format = '{:,.4f}'.format
print('Estimated Coefficients:')
print(lmUD.summary2().tables[1])
print(' ')
print('Number of observations:','{:.0f}'.format(lmUD.nobs),
' Error degrees of freedom:','{:.0f}'.format(lmUD.df_resid))
print('Root Mean Squared Error:',format(np.sqrt(lmUD.mse_resid),'.4g'))
print('R-squared:',format(lmUD.rsquared,'.3g'),' Adjusted R-Squared:',
format(lmUD.rsquared_adj,'.3g'))
print('F-statistic vs. constant model:',format(lmUD.fvalue,'.1f'),
' p-value =',format(lmUD.f_pvalue,'.3g'))
# data indskrives
haemeff0=np.array([75.68,79.58,83.89,85.55,88.11,76.49,84.76,85.58,
86.26,91.43,77.38,86.92,88.85,89.92,93.00])
logKonc0=np.log(np.tile(np.array([0.1,0.25,0.50,0.75,1.00]),3))
gr0=np.repeat(["INH1","INH2","INH3"],5)
# datatabel dannes
# (hvis data indlæses med read_csv har man allerede en datatabel)
datahaemmer=pd.DataFrame({'haemeff':haemeff0,'logKonc':logKonc0,'G':gr0})
print("Opstart er gennemført: haemeff, logKonc og G er indskrevet i tabel datahaemmer")
xxxxxxxxxx
# Opstart ovenfor skal være kørt. Data er i datatabellen datahaemmer
# modeller analyseres
lmUD1=ols(data=datahaemmer,formula='haemeff~G*logKonc').fit()
lmUD2Gskaering=ols(data=datahaemmer,formula='haemeff~G+logKonc').fit()
lmUD3=ols(data=datahaemmer,formula='haemeff~logKonc').fit()
# F-tests laves og udskrives
print('Gruppespecifik regression til gruppespecifik skæring:')
print(anova_lm(lmUD2Gskaering,lmUD1))
print(''); print('Gruppespecifik skæring til simpel lineær regression:')
print(anova_lm(lmUD3,lmUD2Gskaering))
# parametertabel udskrives
print(''); print('Parametertabel for M2Gskaering:')
summaryLM(lmUD2Gskaering)
# konfidensinterval for spredning beregnes
df2=lmUD2Gskaering.df_resid; s22=lmUD2Gskaering.mse_resid
sigmaKonfInt=np.sqrt(df2*s22/st.chi2.ppf(np.array([0.975,0.025]),df2))
print('')
print('Konfidensinterval for spredning under gruppespecifik skæring:',
format(sigmaKonfInt[0],'.2f'),format(sigmaKonfInt[1],'.2f'))