R:在使用自构造的交互变量时无法获得 lme{nlme} 以适应
R: cant get a lme{nlme} to fit when using self-constructed interaction variables
我正在尝试使用自建的交互变量来适应 lme。我需要那些用于 post-hoc 分析。
library(nlme)
# construct fake dataset
obsr <- 100
dist <- rep(rnorm(36), times=obsr)
meth <- dist+rnorm(length(dist), mean=0, sd=0.5); rm(dist)
meth <- meth/dist(range(meth)); meth <- meth-min(meth)
main <- data.frame(meth = meth,
cpgl = as.factor(rep(1:36, times=obsr)),
pbid = as.factor(rep(1:obsr, each=36)),
agem = rep(rnorm(obsr, mean=30, sd=10), each=36),
trma = as.factor(rep(sample(c(TRUE, FALSE), size=obsr, replace=TRUE), each=36)),
depr = as.factor(rep(sample(c(TRUE, FALSE), size=obsr, replace=TRUE), each=36)))
# check if all factor combinations are present
# TRUE for my real dataset; Naturally TRUE for the fake dataset
with(main, all(table(depr, trma, cpgl) >= 1))
# construct interaction variables
main$depr_trma <- interaction(main$depr, main$trma, sep=":", drop=TRUE)
main$depr_cpgl <- interaction(main$depr, main$cpgl, sep=":", drop=TRUE)
main$trma_cpgl <- interaction(main$trma, main$cpgl, sep=":", drop=TRUE)
main$depr_trma_cpgl <- interaction(main$depr, main$trma, main$cpgl, sep=":", drop=TRUE)
# model WITHOUT preconstructed interaction variables
form1 <- list(fixd = meth ~ agem + depr + trma + depr*trma + cpgl +
depr*cpgl +trma*cpgl + depr*trma*cpgl,
rndm = ~ 1 | pbid,
corr = ~ cpgl | pbid)
modl1 <- nlme::lme(fixed=form1[["fixd"]],
random=form1[["rndm"]],
correlation=corCompSymm(form=form1[["corr"]]),
data=main)
# model WITH preconstructed interaction variables
form2 <- list(fixd = meth ~ agem + depr + trma + depr_trma + cpgl +
depr_cpgl + trma_cpgl + depr_trma_cpgl,
rndm = ~ 1 | pbid,
corr = ~ cpgl | pbid)
modl2 <- nlme::lme(fixed=form2[["fixd"]],
random=form2[["rndm"]],
correlation=corCompSymm(form=form2[["corr"]]),
data=main)
第一个模型没有任何问题,而第二个模型给我以下错误:
Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1
到目前为止,我发现的有关此错误的任何信息都没有帮助我解决问题。然而,解决方案可能非常简单。
有人可以帮助我吗?提前致谢!
编辑 1:
当我运行:
modl3 <- lm(form1[["fixd"]], data=main)
modl4 <- lm(form2[["fixd"]], data=main)
总结表明,与 modl3 相比,modl4(具有自建的交互变量)显示出更多的预测变量。所有在 4 中但不在 3 中的都将 NA 显示为系数。因此,问题肯定在于我创建交互变量的方式...
编辑 2:
与此同时,我创建了交互变量 "by hand"(主要是 paste()
和 grepl()
)——它现在似乎可以工作了。但是,我仍然对如何使用 interaction()
函数实现它感兴趣。
我应该只构建最大的交互变量(结合所有 3 个简单变量)。
如果我这样做,模型就会适合。然后可能性彼此非常接近并且系数的数量完全匹配。
我正在尝试使用自建的交互变量来适应 lme。我需要那些用于 post-hoc 分析。
library(nlme)
# construct fake dataset
obsr <- 100
dist <- rep(rnorm(36), times=obsr)
meth <- dist+rnorm(length(dist), mean=0, sd=0.5); rm(dist)
meth <- meth/dist(range(meth)); meth <- meth-min(meth)
main <- data.frame(meth = meth,
cpgl = as.factor(rep(1:36, times=obsr)),
pbid = as.factor(rep(1:obsr, each=36)),
agem = rep(rnorm(obsr, mean=30, sd=10), each=36),
trma = as.factor(rep(sample(c(TRUE, FALSE), size=obsr, replace=TRUE), each=36)),
depr = as.factor(rep(sample(c(TRUE, FALSE), size=obsr, replace=TRUE), each=36)))
# check if all factor combinations are present
# TRUE for my real dataset; Naturally TRUE for the fake dataset
with(main, all(table(depr, trma, cpgl) >= 1))
# construct interaction variables
main$depr_trma <- interaction(main$depr, main$trma, sep=":", drop=TRUE)
main$depr_cpgl <- interaction(main$depr, main$cpgl, sep=":", drop=TRUE)
main$trma_cpgl <- interaction(main$trma, main$cpgl, sep=":", drop=TRUE)
main$depr_trma_cpgl <- interaction(main$depr, main$trma, main$cpgl, sep=":", drop=TRUE)
# model WITHOUT preconstructed interaction variables
form1 <- list(fixd = meth ~ agem + depr + trma + depr*trma + cpgl +
depr*cpgl +trma*cpgl + depr*trma*cpgl,
rndm = ~ 1 | pbid,
corr = ~ cpgl | pbid)
modl1 <- nlme::lme(fixed=form1[["fixd"]],
random=form1[["rndm"]],
correlation=corCompSymm(form=form1[["corr"]]),
data=main)
# model WITH preconstructed interaction variables
form2 <- list(fixd = meth ~ agem + depr + trma + depr_trma + cpgl +
depr_cpgl + trma_cpgl + depr_trma_cpgl,
rndm = ~ 1 | pbid,
corr = ~ cpgl | pbid)
modl2 <- nlme::lme(fixed=form2[["fixd"]],
random=form2[["rndm"]],
correlation=corCompSymm(form=form2[["corr"]]),
data=main)
第一个模型没有任何问题,而第二个模型给我以下错误:
Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1
到目前为止,我发现的有关此错误的任何信息都没有帮助我解决问题。然而,解决方案可能非常简单。
有人可以帮助我吗?提前致谢!
编辑 1:
当我运行:
modl3 <- lm(form1[["fixd"]], data=main)
modl4 <- lm(form2[["fixd"]], data=main)
总结表明,与 modl3 相比,modl4(具有自建的交互变量)显示出更多的预测变量。所有在 4 中但不在 3 中的都将 NA 显示为系数。因此,问题肯定在于我创建交互变量的方式...
编辑 2:
与此同时,我创建了交互变量 "by hand"(主要是 paste()
和 grepl()
)——它现在似乎可以工作了。但是,我仍然对如何使用 interaction()
函数实现它感兴趣。
我应该只构建最大的交互变量(结合所有 3 个简单变量)。
如果我这样做,模型就会适合。然后可能性彼此非常接近并且系数的数量完全匹配。