当将完全随机的两级空模型与包R2MLwiN
拟合时,我收到一条错误消息。我的数据框是联合国儿童基金会为莫桑比克开发的多指标类集调查的子集。我的响应变量agem.18
一个二元("Yes", "No")
,表示该女性是否在 18 岁之前结婚。
> table(moz.20$agem.18, useNA = "ifany")
No Yes
5934 5405
我的两个级别按降序排列为province
和w.id
。这是我运行的代码
# Random effect
F1 <- logit(agem.18) ~ 1 + (1 | province) + (1 | w.id)
rand.eff <- runMLwiN(Formula = F1, D = "Binomial", data = moz.20)
这是我收到的错误消息
Invalid link function specified: NA
Error in runMLwiN(Formula = F1, D = "Binomial", data = df) :
我最初认为logit
函数被car
包掩盖了,但事实并非如此。我也认为问题在链接函数的denom
调用中,但我读到 R2MLwiN 应该自动将denom
创建为一组 1。如果我使用具有相同数据、变量和级别的包lme4
,我不会收到任何错误:
(fit <- glmer(agem.18 ~ 1 + (1 | province), family = binomial("logit"), data = moz.20)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: agem.18 ~ 1 + (1 | province)
Data: moz.20
AIC BIC logLik deviance df.resid
14907.498 14922.170 -7451.749 14903.498 11337
Random effects:
Groups Name Std.Dev.
province (Intercept) 0.5366
Number of obs: 11339, groups: province, 11
Fixed Effects:
(Intercept)
-0.04992
如果我使用包含在包R2MLwiN
demo UserGuide09.R
中的非常相似的公式,则不会遇到同样的问题。
(mymodel1 <- runMLwiN(logit(use) ~ 1 + lc, D = "Binomial", data = bang))
我目前唯一的猜测是,由于某些原因,R2MLwiN
无法将我的响应变量agem.18
识别为二进制。
有什么建议吗?
谢谢
马诺洛
问题解决了。R2MLwiN中似乎有一个错误,正如Chris Charlton在R2MLwiN论坛中指出的那样(见此处的帖子(。我复制了克里斯编写的显示错误的代码。
> library(R2MLwiN)
> data(bang, package = "R2MLwiN")
> (mymodel1 <- runMLwiN(logit(use) ~ 1 + lc, D = "Binomial", data = bang))
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.02) multilevel model (Binomial)
Estimation algorithm: IGLS MQL1 Elapsed time : 0.6s
Number of obs: 2867 (from total 2867) The model converged after 4 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc
Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -1.12288 0.08348 -13.45 3.05e-41 *** -1.28650 -0.95926
lcOne_child 0.93275 0.12156 7.67 1.676e-14 *** 0.69450 1.17100
lcTwo_children 1.09250 0.12509 8.73 2.466e-18 *** 0.84733 1.33768
lcThree_plus 0.87224 0.10302 8.47 2.523e-17 *** 0.67033 1.07416
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
> bang$use.test <- bang$use
> (mymodel1 <- runMLwiN(logit(use.test) ~ 1 + lc, D = "Binomial", data = bang))
Invalid link function specified: NA
Error in runMLwiN(logit(use.test) ~ 1 + lc, D = "Binomial", data = bang) :