Как предсказать игру модели со случайным эффектом в R?

Question

Я работаю над прогнозированием gam модели со случайным эффектом для создания трехмерного графика поверхности по plot_ly.

Вот мой код;

x <- runif(100)
y <- runif(100)
z <- x^2 + y + rnorm(100)
r <- rep(1,times=100) # random effect
r[51:100] <- 2 # replace 1 into 2, making two groups
df <- data.frame(x, y, z, r)

gam_fit <- gam(z ~ s(x) + s(y) + s(r,bs="re"), data = df) # fit

#create matrix data for `add_surface` function in `plot_ly`
newx <- seq(0, 1, len=20)
newy <- seq(0, 1, len=30)
newxy <- expand.grid(x = newx, y = newy)
z <- matrix(predict(gam_fit, newdata = newxy), 20, 30) # predict data as matrix

Однако последняя строка приводит к ошибке;

Error in model.frame.default(ff, data = newdata, na.action = na.act) : 
   variable lengths differ (found for 'r')
In addition: Warning message:
In predict.gam(gam_fit, newdata = newxy) :
  not all required variables have been supplied in  newdata!

Благодаря предыдущему ответу я уверен, что вышеуказанные коды работают без случайного эффекта, как в здесь .

Как я могу предсказать Гам модели со случайным эффектом?

Answer 1

Если вы хотите, чтобы поверхность была обусловлена ​​случайными эффектами (но не для определенного c уровня случайного эффекта), есть два способа.

Первый - предоставить уровень для случайного эффект, но исключить этот термин из прогнозируемых значений, используя аргумент exclude для predict.gam(). Второй - снова использовать exclude, но на этот раз, чтобы не предоставлять никаких данных для случайного эффекта, и вместо этого predict.gam() остановить проверку newdata с помощью аргумента newdata.guaranteed = TRUE.

Опция 1:

newxy1 <- with(df, expand.grid(x = newx, y = newy, r = 2))
z1 <- predict(gam_fit, newdata = newxy1, exclude = 's(r)')
z1 <- matrix(z1, 20, 30)

Вариант 2:

z2 <- predict(gam_fit, newdata = newxy, exclude = 's(r)',
              newdata.guaranteed=TRUE)
z2 <- matrix(z2, 20, 30)

Они дают тот же результат:

> all.equal(z1, z2)
[1] TRUE

Пара примечаний:

1. То, что вы используете, будет зависеть от того, насколько сложна ваша модель. Я бы обычно использовал первый вариант, так как он дает дополнительную проверку на то, что я делаю глупости при создании данных. Но в этом случае, с простой моделью и набором ковариат, кажется достаточно безопасным полагать, что newdata в порядке.

2. В вашем примере используется случайный наклон (это было задумано?) не случайный перехват, так как r не является фактором. Если ваш реальный пример использует факторный случайный эффект, вам нужно быть немного более осторожным при создании newdata, так как вам нужно правильно получить levels коэффициента. Например:

   expand.grid(x = newx, y = newy,
               r = with(df, factor(2, levels = levels(r))))
   

   должен получить правильную настройку для коэффициента r