无法在带有因子变量的 svyglm 中指定起始 glm 参数,R 调查包

Cannot specify starting glm parameters in svyglm with factor variables, R survey package

我正在使用加权分析并使用 svyglm 分析来自复杂加权方案的无响应数据。我想通过指定 binomial(link=log) 作为家庭来拟合一个对数二项式模型来估计适合大多数情况的流行率。但是,在默认拟合器无法找到起始系数集的情况下,我发现在大多数情况下都可以使用的方便的设置是设置 Start <- c(log(mean(response.var)), rep(0, ncov))

当我向 survey 包中的 svyglm 函数提供 start 时,R 抛出一个我似乎无法解析的错误。似乎只要其中一个协变量是一个因素。

示例:

library(survey)
data(api)
apistrat$qmeal <- with(apistrat, cut(meals, quantile(meals)))
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)

还有一个有问题的 GLM 的例子,对一些荒谬的东西进行建模以重现错误:

> svyglm(awards ~ qmeal +emer, family=quasibinomial(link=log), design=dstrat)
Error: no valid set of coefficients has been found: please supply starting values

好的...所以我指定:Start <- c(log(mean(api$awards, na.rm=T)), 0, 0, 0, 0)

> svyglm(awards ~ cut(meals, quantile(meals)) +emer, family=quasibinomial, design=dstrat, start=start)

 > svyglm(awards ~ qmeal +emer, family=quasibinomial(link=log), design=dstrat, start=start)
Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  : 
  length of 'start' should equal 5 and correspond to initial coefs for c("(Intercept)", "qmeal(20.8,39.5]", "qmeal(39.5,69]", "qmeal(69,100]", , "emer")

有趣的是,start 的长度是 5。我进一步注意到 svyglm 始终产生一个额外的 ,(在最后一个 qmeal 变量和 "emer" 之间查看)缺少条目。提供给标准 glm:

时没有这样的问题
glm(awards ~ qmeal +emer, family=quasibinomial(link=log), data=apistrat, start=start)

产生正确的输出:

Call:  glm(formula = awards ~ qmeal + emer, family = quasibinomial(link = log), 
    data = apistrat, start = start)

Coefficients:
     (Intercept)  qmeal(20.8,39.5]    qmeal(39.5,69]     qmeal(69,100]              emer  
        -0.59276           0.13058           0.31311           0.24698          -0.01389  

Degrees of Freedom: 198 Total (i.e. Null);  194 Residual
  (1 observation deleted due to missingness)
Null Deviance:      272.7 
Residual Deviance: 265.7    AIC: NA

svyglm 中调用 glm 的方式似乎有问题。将矢量名称 start 替换为与 svyglm 的参数名称不匹配的任何内容(例如 x)即可解决问题。

您的代码实际上并不 运行,但如果我这样做 start <- c(log(mean(apistrat$awards=="Yes", na.rm=T)), 0, 0, 0, 0)

我确实收到了您询问的错误。发生这种情况是因为对 glm 的调用在设计对象中(有意地)查找其参数,然后在 svyglm 中(并非有意地)查找。正式参数 start 在那里不可见。但是 startstats 包中的一个函数的名称,并且由于复杂的原因[1] 这就是你得到的参数。它的长度不是 5。杂散的逗号是虚假的[2]

一种解决方法是明确指定 start 参数而不是作为变量,因此不需要查找 svyglm(awards ~ qmeal +emer, family=quasibinomial(link=log), design=dstrat, start=c(log(mean(apistrat$awards=="Yes", na.rm=T)), 0, 0, 0, 0))

另一种解决方法是指定一些不是现有函数名称的内容,以便查找进入下一个级别并找到您的变量。例如,这些都对我有用: initial <- c(log(mean(apistrat$awards=="Yes", na.rm=T)), 0, 0, 0, 0) svyglm(awards ~ qmeal +emer, family=quasibinomial(link=log), design=dstrat, start=initial) rose <- c(log(mean(apistrat$awards=="Yes", na.rm=T)), 0, 0, 0, 0) svyglm(awards ~ qmeal +emer, family=quasibinomial(link=log), design=dstrat, start=rose) 我会尝试在下一个版本中修复此问题。

[1] 不,如果没有 运行 一堆实验,我不可能更精确。它们很复杂。

[2] glm.fit 使用 deparse 将名称向量转换为字符串,并且该向量足够长,可以跨两行,这就是逗号的来源。如果变量被称为 m 而不是 qmeal 你就不会得到逗号。您可能会争辩说这是 glm.fit 中的一个错误,但您可能付出的努力超出了它的价值。