Получение стандартизированных коэффициентов для модели glmer?

#r #lme4

#r #lme4

Вопрос:

Меня попросили предоставить стандартизированные коэффициенты для glmer модели, но я не уверен, как их получить. К сожалению, beta функция не работает на glmer моделях:

 Error in UseMethod("beta") : 
  no applicable method for 'beta' applied to an object of class "c('glmerMod', 'merMod')"
  

Есть ли другие функции, которые я мог бы использовать, или мне пришлось бы написать один сам?

Другая проблема заключается в том, что модель содержит несколько непрерывных предикторов (которые работают в аналогичных масштабах) и 2 категориальных предиктора (один с 4 уровнями, один с шестью уровнями). Целью использования стандартизированных коэффициентов было бы сравнить влияние категориальных предикторов на влияние непрерывных, и я не уверен, что стандартизированные коэффициенты являются подходящим способом сделать это. Являются ли стандартизированные коэффициенты приемлемым подходом?

Модель выглядит следующим образом:

 model=glmer(cbind(nr_corr,maximum-nr_corr) ~ (condition|SUBJECT)   categorical_1   categorical_2   continuous_1   continuous_2   continuous_3   continuous_4   categorical_1:categorical_2   categorical_1:continuous_3, data, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)), family = binomial)
  

Ответ №1:

reghelper::beta просто стандартизирует числовые переменные в нашем наборе данных. Итак, предполагая, что ваши катагорические переменные — это factor s, а не числовые фиктивные переменные или другие контрастные кодировки, мы можем довольно просто стандартизировать числовые переменные в нашем наборе данных

 vars <- grep('^continuous(.*)?', all.vars(formula(model)))
f <- function(var, data)
   scale(data[[var]])
data[, vars] <- lapply(vars, f, data = data)
update(model, data = data)
  

Теперь для более общего случая мы можем более или менее просто создать нашу собственную beta.merMod функцию. Однако нам нужно будет принять во внимание, имеет ли смысл стандартизировать y . Например, если у нас есть poisson модель, имеют смысл только положительные целочисленные значения. Кроме того, возникает вопрос, следует ли масштабировать эффекты случайного наклона или нет, и имеет ли смысл задавать этот вопрос в первую очередь. В нем я предполагаю, что категориальные переменные кодируются как character или factor , а не numeric или integer .

 beta.merMod <- function(model, 
                        x = TRUE, 
                        y = !family(model) %in% c('binomial', 'poisson'), 
                        ran_eff = FALSE, 
                        skip = NULL, 
                        ...){
  # Extract all names from the model formula
  vars <- all.vars(form <- formula(model))
  lhs <- all.vars(form[[2]])
  # Get random effects from the 
  ranef <- names(ranef(model))
  # Remove ranef and lhs from vars
  rhs <- vars[!vars %in% c(lhs, ranef)]
  # extract the data used for the model
  env <- environment(form)
  call <- getCall(model)
  data <- get(dname <- as.character(call$data), envir = env)
  # standardize the dataset
  vars <- character()
  if(isTRUE(x))
    vars <- c(vars, rhs)
  if(isTRUE(y))
    vars <- c(vars, lhs)
  if(isTRUE(ran_eff))
    vars <- c(vars, ranef)
  data[, vars] <- lapply(vars, function(var){
    if(is.numeric(data[[var]]))
      data[[var]] <- scale(data[[var]])
    data[[var]]
  })
  # Update the model and change the data into the new data.
  update(model, data = data)
}
  

Функция работает как для линейных, так и для обобщенных линейных моделей со смешанным эффектом (не тестировалась для нелинейных моделей) и используется так же, как и другие бета-функции из reghelper

 library(reghelper)
library(lme4)
# Linear mixed effect model
fm1 <- lmer(Reaction ~ Days   (Days | Subject), sleepstudy)
fm2 <- beta(fm1)
fixef(fm1) - fixef(fm2)
(Intercept)        Days 
  -47.10279   -19.68157 

# Generalized mixed effect model
data(cbpp)
# create numeric variable correlated with period
cbpp$nv <- 
  rnorm(nrow(cbpp), mean = as.numeric(levels(cbpp$period))[as.numeric(cbpp$period)])
gm1 <- glmer(cbind(incidence, size - incidence) ~ nv   (1 | herd),
              family = binomial, data = cbpp)
gm2 <- beta(gm1)
fixef(gm1) - fixef(gm2)
(Intercept)          nv 
  0.5946322   0.1401114
  

Однако обратите внимание, что в отличие beta от функции возвращает обновленную модель, а не сводную информацию о модели.

Другая проблема заключается в том, что модель содержит несколько непрерывных предикторов (которые работают в аналогичных масштабах) и 2 категориальных предиктора (один с 4 уровнями, один с шестью уровнями). Целью использования стандартизированных коэффициентов было бы сравнить влияние категориальных предикторов на влияние непрерывных, и я не уверен, что стандартизированные коэффициенты являются подходящим способом сделать это. Являются ли стандартизированные коэффициенты приемлемым подходом?

Теперь это отличный вопрос, и он лучше подходит для stats.stackexchange него, но я не уверен в ответе на него.

Комментарии:

1. Большое тебе спасибо, Оливер! Я все еще пытаюсь разобраться во всем этом. Могу ли я спросить, что именно делает строка ‘cbpp $ nv <- rnorm(nrow (cbpp), mean = as.numeric(levels(cbpp $ period)) [as.numeric (cbpp $ period)])’? Насколько я понимаю, он генерирует новую переменную, состоящую из случайных чисел с определенным средним значением, но я не уверен, что именно делает команда для среднего значения.

2. Вы правильно поняли. Данные примера, которые я использовал cbpp , не имеют числовой переменной, и у меня не было примера с одной под рукой. Но у него есть коэффициент period (значения: от «1» до «4»). period используется как фиксированный эффект в примерах ?glmer , поэтому я создал числовой прокси для фактора, сгенерировав некоторые случайные данные со средним значением от 1 до 4 (эквивалент периода). Короче говоря, это rnorm(nrow(cbpp), mean = (1:4)[match with period as numeric]) говорит. rnorm(n, mean = mu) просто генерирует случайные n числа со средними значениями mu ( mu может быть больше 1 числа или вектора длины n)

3. Я мог бы также использовать cbpp$nv <- rnorm(nrow(cbpp), mean = as.numeric(cbpp$period)) .

Ответ №2:

Еще раз, большое спасибо, Оливер! Для всех, кого интересует ответ на последнюю часть моего вопроса,

Другая проблема заключается в том, что модель содержит несколько непрерывных предикторов (которые работают в аналогичных масштабах) и 2 категориальных предиктора (один с 4 уровнями, один с шестью уровнями). Целью использования стандартизированных коэффициентов было бы сравнить влияние категориальных предикторов на влияние непрерывных, и я не уверен, что стандартизированные коэффициенты являются подходящим способом сделать это. Являются ли стандартизированные коэффициенты приемлемым подходом?

вы можете найти ответ здесь. Проблема в том, что использование стандартизированных коэффициентов регрессии в любом случае не лучший подход для смешанных моделей, не говоря уже о таком, как у меня…