时间戳

首页 > 解决方案 > 将 glmer 输出（logit 回归）转换为概率

r - 将 glmer 输出（logit 回归）转换为概率

问题描述

我有这样的数据。我在所有 Q 变量上运行 glm 。

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
")

这是我用来提取我想要的系数的循环。

# 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

这个循环工作正常，但我需要将输出（系数）从对数赔率转换为概率。关于如何做到这一点的任何建议？

标签： rregression

解决方案

坏消息：实际上没有任何明智的方法可以将逻辑回归的系数（在对数优势比或 logit 尺度上）转换为概率尺度。从对数赔率到概率的转换取决于基线水平，因此要获得概率，您必须对特定情况下的概率进行预测：参见例如这个 CrossValidated question。

好消息：对系数取幂可以得到优势比，这通常比对数优势比更容易理解并且可以说更容易理解。

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 值等的估计值。还有其他选项，例如conf.int=TRUE您是否需要表格中的置信区间。您可以使用任何 tidyverse 工具来操作它（例如%>% filter(term=="condition")，如果您对截取不感兴趣）。

这个例子中的许多答案都是假的，但那是因为你的数据集太小了......我希望你的真实数据集比这个大......

解释为什么您通常不能将优势比转换为概率（不指定基线）实际上更像是一个统计/交叉验证问题，但我将根据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()是用于逆 logit（逻辑）变换的内置 R 函数。

转换截距参数确实有意义：它给出了具有基线特征的个体的预测概率；由于我们已将预测变量居中，因此这对应于具有总体平均 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 的概率变化会有所不同。

如果需要，您可以对 GRE 系数进行逻辑转换：

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

不过，这意味着什么？如果您随后将其 GRE 提高 1 个标准差，则它是基线对数几率为零的个人（不是具有平均 GPA 和 GRE 的个人）的成功概率。不是一件容易解释的事情......

还有其他经验法则/近似规则可以理解对数赔率的含义，例如除以 4 规则，但它们都以某种方式依赖于指定基线水平。

推荐阅读