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