Я пытался вписать регрессионную модель штрафованного сплайна в R, используя связь между штрафованными сплайнами и линейными смешанными моделями. Хотя я достаточно знаком с функцией R lme, использование ее конкурента lmer представляет некоторые трудности. Вот игрушечный пример, который я хотел бы обобщить до lmer:
require(nlme)
grid <- seq(0, 1, len = 100)
y <- rep(0, length(grid))
for(i in 1:length(grid)){
y[i] <- sin(3*pi*grid[i]) + rnorm(1, 0, 1)
}
X <- cbind(rep(1, length(grid)), grid, grid^2, grid^3)
knots <- seq(0.01, 0.99, len = 100)
Z <- outer(grid, knots, "-")
Z <- Z*(Z>0)
Z <- Z^3
data.mixed <- data.frame(X, Z)
data.mixed$all <- rep(1 ,nrow(data.mixed))
fit.mixed <- lme(y~ X-1, random = list( all = pdIdent(~Z-1) ), data = data.mixed )
curve(sin(3*pi*x), 0, 1, lwd = 2, ylim = c(-1.5, 1.5))
beta.hat <- fit.mixed$coefficients$fixed
u.hat <- unlist(fit.mixed$coefficients$random)
f.hat <- X%*%beta.hat + Z%*%u.hat
lines(grid, f.hat, lwd = 2, col = "red")
Сложная часть - это функция pdIdent, которая указывает, что, хотя у нас есть несколько случайных эффектов, их ковариационная структура кратна единичной матрице; lme затем оценивает это скалярное кратное.
Если кто-нибудь знаком с lmer, я был бы очень признателен за некоторые советы о том, как выполнить это с помощью этой функции.
Заранее спасибо