2

我有两个没有随机效应的回归模型：一个是使用 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)")

有什么建议么？

4

1 回答 1

2

试试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）。并且很容易为此添加随机拦截。

于 2015-04-14T15:50:49.457 回答