问题描述
我的任务是像下面的SAS代码那样在python中计算相同的最低有效差（LSD）。我对统计数据不太熟悉，因此不胜感激任何公式或代码建议。
proc mixed covtest CL;
by loc;
class line rep;
model &var=line ;
random rep;
estimate 'testmean' intercept 1;
lsmean line/diff=all;
ods output diffs=db30;
ods output lsmeans;
ods listing exclude all;
proc means data=all;
by loc;
var &var;
output out=minmax max=max min=min mean=mean std=std;
proc print data=minmax;
var loc _Freq_ mean std min max;
title 'Minimum and Maximum values by location on Raw Data';
这是我在python中尝试过的方法，但是返回的值太低。 new_data是一个具有4列的数据集-一列每个品种的条目代码，一列代表Rep，一列代表Plot，一列包含Yield数据：
# Calculate the lsd
lsd = gs.calclsd(new_data,'YLD',"ENTRYCODE")
def calclsd(dataset,depvar,indepvar1):
""""
#############################################
# Purpose: #
# Calculate the Least Squares Difference #
# of a balanced dataset (no missing values) #
# Inputs: #
# dataSet - test dataset (must be balanced) #
# depVar - dependent variable #
# indepVar1 - independent variable 1 #
#############################################
Get the ordinary least squares
modelLSD = ols(formula,data=dataSet).fit()
calculate the anova
"""
"""aov_table = anova_lm(modelLSD)"""
"""one-way anova"""
model=ols(str(depvar) + ' ~ C(ENTRYCODE)',data=dataset).fit()
print(model.summary())
aov_table=smi.stats.anova_lm(model,typ=3)
aov_table=anova_table(aov_table)
print(aov_table)
degreesf = aov_table.at['C(ENTRYCODE)','df']
meansq = aov_table.at['C(ENTRYCODE)','mean_sq']
"""calculate LSD"""
n = len(dataset['DESIGNATION'].unique())
mean = dataset[depvar].mean()
prob = 0.05
stats.ttest(dataset[depvar])
print(tval)
print(pval)
#tval = t.ppf(prob,degreesf)
lsd = round(tval * math.sqrt(meansq * (2 / n)),1)
return lsd
def anova_table(anov):
anov['mean_sq'] = anov[:]['sum_sq']/anov[:]['df']
anov['eta_sq'] = anov[:-1]['sum_sq']/sum(anov['sum_sq'])
anov['omega_sq'] = (anov[:-1]['sum_sq']-(anov[:-1]['df']*anov['mean_sq'][-1]))/(sum(anov['sum_sq'])+anov['mean_sq'][-1])
cols = ['sum_sq','df','mean_sq','F','PR(>F)','eta_sq','omega_sq']
anov = anov[cols]
return anov
再次，任何帮助都是很好的，因为我在这个问题上花了很多时间。
干杯！
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）