r - R 中系数相乘的非线性随机效应回归

Question

我有两个没有随机效应的回归模型:一个是使用 lm 的 OLS，另一个包括使用 nle 的系数乘法。 我希望为两者添加个人级别的随机效应。我已经设法使用 lme4 包对 OLS 函数执行此操作，但无法找到对乘法模型执行此操作的方法。

以下代码生成一个与我正在处理的结构类似的数据集:

df <- data.frame(id = rep(1:1000, each=10), jit = rep(rnorm(1000, 0, 0.2), each = 10), a = sample(1:5, 10000, T), b = sample(1:5, 10000,T), c = sample(1:5, 10000, T))
df <- cbind(df, model.matrix(~ as.factor(a) + as.factor(b) + as.factor(c), data.frame(rbind(as.matrix(df), t(matrix(rep(1:5, each = 5), nrow=5)))))[1:nrow(df),2:13])
colnames(df)[6:17] <- (dim_dummies <- as.vector(outer(2:5, letters[1:3], function(x, y) paste(y, x, sep=""))))
true_vals <- list(vL2 = 0.4, vL3 = 0.5, vL4 = 0.8, vA = 0.7, vB = 1.1, vC = 0.9)
attach(df)
attach(true_vals)
df$val <- 
  (a2 * vA + b2*vB + c2*vC) * vL2 + 
  (a3 * vA + b3*vB + c3*vC) * vL3 + 
  (a4 * vA + b4*vB + c4*vC) * vL4 + 
  (a5 * vA + b5*vB + c5*vC) + runif(1, -.2, .2) + jit
detach(true_vals)
detach(df)

df[1:15, ]
   id      jit a b c a2 a3 a4 a5 b2 b3 b4 b5 c2 c3 c4 c5     val
1   1 -0.14295 4 4 1  0  0  1  0  0  0  1  0  0  0  0  0  1.1698
2   1 -0.14295 5 1 4  0  0  0  1  0  0  0  0  0  0  1  0  1.1498
3   1 -0.14295 5 4 4  0  0  0  1  0  0  1  0  0  0  1  0  2.0298
4   1 -0.14295 5 1 5  0  0  0  1  0  0  0  0  0  0  0  1  1.3298
5   1 -0.14295 5 4 2  0  0  0  1  0  0  1  0  1  0  0  0  1.6698
6   1 -0.14295 1 5 1  0  0  0  0  0  0  0  1  0  0  0  0  0.8298
7   1 -0.14295 3 2 5  0  1  0  0  1  0  0  0  0  0  0  1  1.4198
8   1 -0.14295 3 2 1  0  1  0  0  1  0  0  0  0  0  0  0  0.5198
9   1 -0.14295 3 2 4  0  1  0  0  1  0  0  0  0  0  1  0  1.2398
10  1 -0.14295 5 3 3  0  0  0  1  0  1  0  0  0  1  0  0  1.4298
11  2 -0.01851 4 5 3  0  0  1  0  0  0  0  1  0  1  0  0  1.9643
12  2 -0.01851 2 1 3  1  0  0  0  0  0  0  0  0  1  0  0  0.5843
13  2 -0.01851 2 1 3  1  0  0  0  0  0  0  0  0  1  0  0  0.5843
14  2 -0.01851 1 1 1  0  0  0  0  0  0  0  0  0  0  0  0 -0.1457
15  2 -0.01851 2 3 1  1  0  0  0  0  1  0  0  0  0  0  0  0.6843

...

a、b 和 c 代表三个 1:5 维度尺度上的分数。 a2 到 c5 是虚拟变量，代表相同尺度的 2:5 水平。每个个体 (id) 有 10 个观察值。 val 是我希望使用回归模型预测的分数的代理。 (但是，实际数据中的值可能与此处的结构不对应。)

我有两个没有随机效应的回归模型。一种是使用 12 个虚拟变量作为 val 预测变量的常规 OLS:

additive.formula <- as.formula("val ~ 
  a2 + a3 + a4 + a5 + 
  b2 + b3 + b4 + b5 + 
  c2 + c3 + c4 + c5")
fit.additive <- lm(additive.formula, data = df)

第二个假设级别之间的相对距离对于三个维度 (a,b,c) 是共享的，但维度在比例方面有所不同。剩下 6 个系数(cA、cB、cC、cL2、cL3、cL4)+ 截距。

multiplicative.formula <- as.formula(" val ~ intercept +
  (a2 * cA + b2*cB + c2*cC) * cL2 + 
  (a3 * cA + b3*cB + c3*cC) * cL3 + 
  (a4 * cA + b4*cB + c4*cC) * cL4 + 
  (a5 * cA + b5*cB + c5*cC)")
multiplicative.start <- list(intercept = 0, cA = 1, cB = 1, cC = 1, cL2 = 1, cL3 = 1, cL4 = 1)
fit.multiplicative <- nls(multiplicative.formula, start=multiplicative.start, data=df, control = list(maxiter = 5000))

由于每个人有 10 个观察值，我们不能期望它们完全独立。因此，我希望在变量 id 定义的个体级别添加随机效应。我找到了一种使用 lme4 包来做到这一点的方法:

require(lme4)
additive.formula.re <- as.formula("val ~ (1 | id) +
  a2 + a3 + a4 + a5 + 
  b2 + b3 + b4 + b5 + 
  c2 + c3 + c4 + c5")
fit.additive.re <- lmer(additive.formula.re, data=df)

问题是是否可以使用类似于乘法模型的回归模型(也许使用 lme4 或 nlme 包)对 id 变量添加随机效应？公式应该类似于

multiplicative.formula.re <- as.formula(" val ~ (1 | id) + intercept +
  (a2 * cA + b2*cB + c2*cC) * cL2 + 
  (a3 * cA + b3*cB + c3*cC) * cL3 + 
  (a4 * cA + b4*cB + c4*cC) * cL4 + 
  (a5 * cA + b5*cB + c5*cC)")

有什么建议吗？

Answer 1

尝试nlme。这应该是你所需要的(如果我理解正确的话):

library(nlme)
fit.multiplicative.nlme <- nlme( model = val ~ intercept +
                                   (a2 * cA + b2*cB + c2*cC) * cL2 + 
                                   (a3 * cA + b3*cB + c3*cC) * cL3 + 
                                   (a4 * cA + b4*cB + c4*cC) * cL4 + 
                                   (a5 * cA + b5*cB + c5*cC),
                                 fixed = intercept + cA +cB + cC + cL2 + cL3 + cL4 ~ 1,
                                 random = intercept ~ 1|id,
                                 start = unlist(multiplicative.start), data=df)

但是，当我使用您提供的不可重现的数据进行尝试时，这并没有收敛(您应该设置随机种子)。您可以在nlmeControl中尝试不同的设置。

---

以下内容不正确:

我没有看到非线性最小二乘的原因。让我们恢复虚拟编码:

df$id1 <- seq_len(nrow(df))
df$a1 <- as.integer(rowSums(df[, paste0("a", 2:5)]) == 0)
df$b1 <- as.integer(rowSums(df[, paste0("b", 2:5)]) == 0)
df$c1 <- as.integer(rowSums(df[, paste0("c", 2:5)]) == 0)
library(reshape2)
DFm <- melt(df, id.vars = c("id", "jit", "a", "b", "c", "val", "id1"))
DFm <- DFm[DFm$value == 1,]
DFm$g <- paste0("fac", substr(DFm$variable, 1, 1))
DF <- dcast(DFm, ... ~ g, value.var = "variable")


fit1 <- lm(val ~ faca + facb + facc, data = DF)

#compare results:
coef(fit.multiplicative)
prod(coef(fit.multiplicative)[c("cA", "cL2")])
coef(fit1)["facaa2"]
prod(coef(fit.multiplicative)[c("cA", "cL3")])
coef(fit1)["facaa3"]

如您所见，这基本上是相同的模型(差异是由于 nls 内的数值优化造成的)。并且很容易为此添加随机截距。