Использование пакета эффектов внутри функции

Question

У меня большой кадр с множеством переменных, которые я собираюсь проанализировать таким же образом. В частности, я хочу построить доверительные интервалы эффекта в модели смешанного эффекта. Я хочу написать функцию, которая делает пользовательский график для одной зависимой переменной. Непосредственное применение функции effect () проходит хорошо. Но тот же код внутри функции вызывает ошибку.

Я попробовал два варианта функции. Оба вызывают ошибки.

Вот мой воспроизводимый пример:

library(nlme)
library(effects)

df <- data.frame(y = rnorm(90), x = gl(3, 30), b = factor(rep(1:30, 3)))

fit <- lme(fixed = y ~ x, random = ~ 1 | b, data = df, method = "REML")
ef <- effect("x", fit)
bp <- barplot(as.vector(ef$fit), col = c("tomato", "skyblue", "limegreen"), 
        ylim = c(min(ef$lower), max(ef$upper) + (max(ef$upper) - min(ef$lower)) * 0.2 ))
arrows(x0 = bp, y0 = ef$lower, y1 = ef$upper, code = 3, angle = 90)

test1 <- function(y, x, b)
{
  fit <- lme(fixed = y ~ x, random = ~ 1 | b, method = "REML")
  ef <- effect("x", fit)
  bp <- barplot(as.vector(ef$fit), col = c("tomato", "skyblue", "limegreen"), 
                ylim = c(min(ef$lower), max(ef$upper) + (max(ef$upper) - min(ef$lower)) * 0.2 ))
  arrows(x0 = bp, y0 = ef$lower, y1 = ef$upper, code = 3, angle = 90)
}

test1(df$y, df$x, df$b)
# Error in eval(predvars, data, env) : object 'y' not found 

test2 <- function(y, x, b)
{
  frame <- data.frame(y, x, b)

  fit <- lme(fixed = y ~ x, random = ~ 1 | b, frame, method = "REML")
  ef <- effect("x", fit)
  bp <- barplot(as.vector(ef$fit), col = c("tomato", "skyblue", "limegreen"), 
                ylim = c(min(ef$lower), max(ef$upper) + (max(ef$upper) - min(ef$lower)) * 0.2 ))
  arrows(x0 = bp, y0 = ef$lower, y1 = ef$upper, code = 3, angle = 90)
}

test2(df$y, df$x, df$b)
# Error in as.data.frame.default(data, optional = TRUE) : 
#  cannot coerce class ‘"function"’ to a data.frame 

Answer 1

Проще:

function(df) {
  fit <- lme(fixed = y ~ x, random = ~ 1 | b, data = df, method = "REML")
  ef <- effect("x", fit)
  bp <- barplot(as.vector(ef$fit), col = c("tomato", "skyblue", "limegreen"),
                ylim = c(min(ef$lower), 
                         max(ef$upper) + (max(ef$upper) - min(ef$lower)) * 0.2 ))
  arrows(x0 = bp, y0 = ef$lower, y1 = ef$upper, code = 3, angle = 90)
}

Вам нужно передать data в lme, формула фактически не передает никаких данных.

Тем не менее, ваш test2 долженработай. Я могу повторить вашу ошибку, но это действительно очень странно. Каким-то образом код работает в глобальной среде, но не в замыкании. Очень удивительно.