使用r中的Survey包进行单向方差分析

Question

我正试图确定在复杂的调查设计上运行单向方差分析的最佳方法． 在仔细阅读了Lumley的Survey包文档后，我一点也不知道。

survey::anova函数的意思是"拟合和比较复杂调查数据的分层对数线性模型"，这不是我正在做的。

我想做什么我收集了一个分类自变量[3个水平]和一个定量因变量的数据。我想用方差分析来检查因变量是否根据自变量的水平而变化。

下面是我的过程的一个例子:

加载调查包并创建复杂的调查设计对象

library(survey)
df <- data.frame(sex = c('F', 'O', NA, 'M', 'M', 'O', 'F', 'F'),
married = c(1,1,1,1,0,0,1,1),
pens = c(0, 1, 1, NA, 1, 1, 0, 0),
weight = c(1.12, 0.55, 1.1, 0.6, 0.23, 0.23, 0.66, 0.67))
svy_design <- svydesign(ids=~1, data=df, weights=~weight)

借用这篇文章，

方法一:使用survey::aov

summary(aov(weight~sex,data = svy_design))

然而，我得到一个错误说:

Error in h(simpleError(msg, call)) : 
error in evaluating the argument 'object' in selecting a method for function 'summary': object 'api00' not foun

方法2:使用调查::glm代替方差分析

同一篇文章有一个反对使用方差分析的答案/解释:

根据我们研究所的主要统计学家的说法，在任何通用的建模环境中都不容易实现这种分析。原因是ANOVA和ANCOVA是线性模型，在70年代通用线性模型(后来的广义线性模型- GLMs)出现后没有得到进一步发展。正常的线性回归模型产生的结果实际上与方差分析相同，但在变量选择方面更加灵活。由于存在glm的加权方法(见R中的调查包)，因此实际上不需要在ANOVA中开发分层抽样设计的加权方法……

summary(svyglm(weight~sex,svy_design))

我得到了这样的输出:

call:
svyglm(formula = weight ~ sex, design = svy_design)
Survey design:
svydesign(ids = ~1, data = df, weights = ~weight)
Coefficients:
Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.8730     0.1478   5.905  0.00412 **
sexM         -0.3756     0.1855  -2.024  0.11292   
sexO         -0.4174     0.1788  -2.334  0.07989 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.04270091)
Number of Fisher Scoring iterations: 2

<标题>我的问题:

1. 为什么方法1抛出错误?
2. 是否可以使用survey::aov函数来完成我的目标?
3. 如果我要使用survey::glm[方法2]，我应该看哪个值来识别平均值的差异?它是截距的p值吗?

我不是一个统计迷，请尽量用最简单的术语解释。谢谢你! !

Answer 1

没有survey::aov这样的函数，所以您不能使用它来完成您的目标。您的代码使用stats::aov

可以使用survey::svyglm。我将使用包中的一个示例，因此我可以实际运行代码

> model<-svyglm(api00~stype, design=dclus2)
> summary(model)
Call:
svyglm(formula = api00 ~ stype, design = dclus2)
Survey design:
dclus2<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)    
(Intercept)   692.81      30.28  22.878  < 2e-16 ***
stypeH        -94.47      27.66  -3.415  0.00156 ** 
stypeM        -50.46      23.01  -2.193  0.03466 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 17528.44)
Number of Fisher Scoring iterations: 2

有三种学校类型，E,M,H。这里的两个系数估计了E的平均值与其他两组的平均值之间的差异，p值检验了H和E具有相同平均值，M和E具有相同平均值的假设。

如果您想对三组之间的均值差异进行总体测试，您可以使用regTermTest函数，该函数测试模型中的一个或一组术语，例如，

> regTermTest(model,~stype)
Wald test for stype
in svyglm(formula = api00 ~ stype, design = dclus2)
F =  12.5997  on  2  and  37  df: p= 6.7095e-05 

该F检验类似于stats::aov给出的F检验。不一样，因为这是调查数据