Нелинейная регрессия случайных эффектов с умножением коэффициентов в R

Question

У меня есть две модели регрессии без случайных эффектов: одна из них - OLS с использованием lm, другая - умножение коэффициентов с использованием nle. Я хочу добавить индивидуальные случайные эффекты для обоих. Я смог сделать это для функции OLS, используя пакет lme4, но не смог найти способ сделать это для мультипликативной модели.

Следующий код создает набор данных с аналогичной структурой в один я работаю:

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 

...

а, Ь, с представляют результаты на три 1: 5 размерных шкал , a2 - c5 - фиктивные переменные, представляющие уровни 2: 5 в одних и тех же масштабах. Есть 10 наблюдений на человека (id). val является прокси для оценки, которую я хочу предсказать с использованием моделей регрессии. (Однако значения в фактических данных могут не соответствовать структуре здесь.)

У меня есть две модели регрессии без случайных эффектов. Одним из них является регулярное МНК с использованием 12 фиктивных переменных, как предикторы Валу:

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

Второй предполагает, что относительное расстояние между уровнями является общим для всех трех dinensions (а, б, в), но, что размеры отличаются по масштабу. Это оставляет 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 наблюдений, мы не можем ожидать, что они будут полностью независимыми. Поэтому я хочу добавить случайный эффект на уровне индивидуума, как определено идентификатором переменной. Я нашел способ сделать это с пакетом 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 пакетов? Формула должна выглядеть примерно так:

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). И легко добавить случайный перехват к этому.