Преобразуйте вывод glmer (lo git регрессия) в вероятности

Question

У меня есть такие данные. Я использую glm для всех переменных Q.

dat <- read_table2("condition   school  Q5_3    Q6  Q7_1    Q7_2    Q7_3    Q7_4    Q13_1   Q13_2   Q13_3
0   A   1   1   1   1   1   1   0   1   1
1   B   1   0   0   NA  NA  NA  NA  1   1
1   C   1   0   1   1   1   1   0   1   1
1   A   0   0   0   NA  NA  NA  NA  1   1
1   B   1   0   0   NA  NA  NA  NA  1   1
0   C   1   1   1   1   1   0   0   0   0
0   A   0   0   0   NA  NA  NA  NA  NA  NA
0   B   1   1   1   1   1   1   1   1   1
0   C   1   1   0   NA  NA  NA  NA  1   0
0   A   1   0   0   NA  NA  NA  NA  1   0
0   B   1   0   1   1   0   1   1   NA  NA
0   C   1   0   1   1   1   1   1   1   0
1   A   1   1   1   1   0   1   0   1   1
1   B   0   0   0   NA  NA  NA  NA  1   1
0   C   1   0   0   NA  NA  NA  NA  NA  NA
")

Это l oop, который я использую для получения нужных мне коэффициентов.

# We only need the condition and school
# Apply
models <- function(x)
{
  model1 <- glmer(x~ (1|school) + condition ,data=dat , family = binomial, na.action = na.exclude)
  return(model1)
}


y <- apply(dat[,-c(1,2)],2,models)
#Extract results
extract <- function(x)
{
  z <- as.data.frame(summary(x)$coefficient)
  z$id <- rownames(z)
  z <- z[,c(dim(z)[2],1:(dim(z)[2]-1))]
  rownames(z)<-NULL
  return(z)
}
#Extract summary with function
DF <- as.data.frame(do.call(rbind,lapply(y,extract)))
#Format variables
DF$var <- gsub("\\..*","",rownames(DF))
#Arrange columns
DF_glm <- DF[,c(dim(DF)[2],1:(dim(DF)[2]-1))]
rownames(DF)<-NULL

Этот цикл работает нормально, но мне нужно преобразовать выходные данные (коэффициенты) из логарифма шансов в вероятности. Есть предложения, как это сделать?

Answer 1

Плохие новости : на самом деле нет никакого разумного способа преобразовать коэффициенты логистической c регрессии (которые находятся в логарифмическом соотношении шансов или шкале lo git) в шкалу вероятностей. Преобразование из логарифмических шансов в вероятности зависит от базового уровня, поэтому для получения вероятностей вам нужно будет сделать прогнозы вероятностей для определенных c случаев: см., Например, этот вопрос с перекрестной проверкой .

Хорошие новости : возведение в степень коэффициентов дает шансы крысы ios, которые, как правило, хорошо понятны и, возможно, легче понять, чем логарифм отношения шансов.

library(broom.mixed)
dd <- dat[,-c(1,2)]
## find (and drop) examples with no variation
uu <- apply(dd,2,function(x) length(unique(na.omit(x))))
modList <- apply(dd[,uu>1],2,models)
## generate list of models
purrr:::map_dfr(modList,tidy,
        effect="fixed",
        exponentiate=TRUE,.id="Q")

Это дает вам таблицу (таблица) с оценками по шкале отношения шансов , стандартными ошибками, значениями p и т. Д. c. Если вам нужны доверительные интервалы в таблице, есть и другие варианты, например conf.int=TRUE. Вы можете управлять им с помощью любого из инструментов tidyverse (например, %>% filter(term=="condition"), если вас не интересуют перехваты).

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

Объяснение того, почему вы не можете обычно преобразовать шансы крысы ios в вероятности (без указания базового уровня), на самом деле является скорее статистическим вопросом / CrossValidated , но я приведу короткий пример, основанный на сайт статистики UCLA

• Импорт данных: масштабируйте переменные-предикторы для GRE и GPA, чтобы получить более интерпретируемые значения параметров.

library(tidyverse)
dd <- (haven::read_dta("https://stats.idre.ucla.edu/stat/stata/dae/binary.dta")
    %>% mutate_at(c("gre","gpa"), ~drop(scale(.)))
)

• Подобрать модель и извлечь коэффициенты

m <- glm(admit~gre+gpa, family=binomial, dd)
cc <- coef(m)
## (Intercept)         gre         gpa 
##  -0.8097503   0.3108184   0.2872088

• преобразование:

plogis() - встроенная функция R для обратного lo git (logisti c) преобразование.

Преобразование параметра перехвата имеет смысл: оно дает прогнозируемую вероятность для человека с базовыми характеристиками; поскольку мы центрировали предикторы, это соответствует индивиду со средним средним GPA и GRE для популяции.

int_prob <- plogis(cc["(Intercept)"])  ## 0.307

Мы также можем предсказать вероятность для человека со средними GRE и GPA на одно стандартное отклонение выше среднее значение (единицы измерения параметра GPA - это «за стандартное отклонение», потому что мы масштабировали переменную GPA по ее стандартному отклонению):

gre_prob <- with(as.list(cc), plogis(`(Intercept)`+gre)) ## 0.3777

Мы могли бы вычислить разницу между этими прогнозами, что является одним из способов определения влияние GRE на шкалу вероятности:

gre_prob-int_prob ## 0.0698

Однако это применимо только для этого конкретного сравнения (человек со средним GPA и GRE на 1 SD выше среднего по сравнению с человеком со средним GPA и GRE ). Изменение вероятности на единицу GRE было бы другим, если бы мы начали с другой базовой линии или , сделав прогноз для другого изменения GRE.

Вы можете логистически c преобразовать коэффициент GRE если вы хотите:

plogis(cc["gre"])  ## 0.577

Но что это значит? Это вероятность успеха для человека с исходными логарифмическими шансами, равными нулю (что составляет не человека со средним GPA и GRE), если вы затем увеличите его GRE на 1 стандартное отклонение. Не то, что легко объяснить ...

Существуют и другие практические правила / аппроксимации для понимания значения log-odds-rat ios, например, правило деления на 4 , но все они так или иначе зависят от определения базового уровня.

Answer 2

Вы можете попробовать это. Вы получите предупреждения из-за ваших данных:

library(lme4)
#Preprocess data
#If you omit NA variables will be constants so that the model
#can no be fitted and will produce error
#I set NA to zero in order to get models working
#Please check your data
dat[is.na(dat)] <- 0
# We only need the condition and school
# Apply
models <- function(x)
{
  model1 <- glmer(x~ (1|school) + condition ,data=dat , family = binomial, na.action = na.exclude)
  return(model1)
}

y <- apply(dat[,-c(1,2)],2,models)
#Extract results and we will extract the logs
extract <- function(x)
{
  z <- as.data.frame(summary(x)$coefficient)
  z$id <- rownames(z)
  z <- z[,c(dim(z)[2],1:(dim(z)[2]-1))]
  z$odds <- exp(z$Estimate)
  z$prob <- z$odds / (1 + z$odds)
  rownames(z)<-NULL
  return(z)
}

#Extract summary with function
DF <- as.data.frame(do.call(rbind,lapply(y,extract)))
#Format variables
DF$var <- gsub("\\..*","",rownames(DF))
#Arrange columns
DF_glm <- DF[,c(dim(DF)[2],1:(dim(DF)[2]-1))]
rownames(DF)<-NULL

Вы получите это:

            id      Estimate   Std. Error       z value   Pr(>|z|)         odds
1  (Intercept)  2.079442e+00 1.060660e+00  1.960517e+00 0.04993534 8.000000e+00
2    condition -1.386294e+00 1.369306e+00 -1.012407e+00 0.31134363 2.500000e-01
3  (Intercept) -2.231436e-01 6.708203e-01 -3.326428e-01 0.73940393 8.000000e-01
4    condition -1.386294e+00 1.284523e+00 -1.079229e+00 0.28048568 2.500000e-01
5  (Intercept)  2.231436e-01 6.708203e-01  3.326428e-01 0.73940394 1.250000e+00
6    condition -9.162907e-01 1.095445e+00 -8.364553e-01 0.40289882 4.000000e-01
7  (Intercept)  2.231436e-01 6.708203e-01  3.326428e-01 0.73940394 1.250000e+00
8    condition -9.162907e-01 1.095445e+00 -8.364553e-01 0.40289882 4.000000e-01
9  (Intercept) -2.231436e-01 6.708204e-01 -3.326428e-01 0.73940397 8.000000e-01
10   condition -1.386294e+00 1.284523e+00 -1.079229e+00 0.28048583 2.500000e-01
11 (Intercept) -2.231436e-01 6.708204e-01 -3.326428e-01 0.73940395 8.000000e-01
12   condition -4.700036e-01 1.095445e+00 -4.290527e-01 0.66788485 6.250000e-01
13 (Intercept) -7.440587e-01 1.454336e+00 -5.116141e-01 0.60892109 4.751814e-01
14   condition -5.938497e+04 2.739708e+07 -2.167566e-03 0.99827053 0.000000e+00
15 (Intercept)  2.231435e-01 6.708204e-01  3.326427e-01 0.73940398 1.250000e+00
16   condition  3.442056e+01 1.351269e+07  2.547276e-06 0.99999797 8.884999e+14
17 (Intercept) -1.252763e+00 8.017837e-01 -1.562470e+00 0.11817732 2.857143e-01
18   condition  3.800559e+01 2.739708e+07  1.387213e-06 0.99999889 3.203452e+16
        prob   var
1  0.8888889  Q5_3
2  0.2000000  Q5_3
3  0.4444444    Q6
4  0.2000000    Q6
5  0.5555556  Q7_1
6  0.2857143  Q7_1
7  0.5555556  Q7_2
8  0.2857143  Q7_2
9  0.4444444  Q7_3
10 0.2000000  Q7_3
11 0.4444444  Q7_4
12 0.3846154  Q7_4
13 0.3221173 Q13_1
14 0.0000000 Q13_1
15 0.5555556 Q13_2
16 1.0000000 Q13_2
17 0.2222222 Q13_3
18 1.0000000 Q13_3