如何从GLM输出中获得概率

Question

我现在非常困难,因为我试图弄清楚如何从glmR中的输出计算概率.我知道数据非常微不足道但我真的很想展示如何从输出中获得概率这个.我正在考虑尝试,inv.logit()但不知道括号内放置了哪些变量.

数据来自占用率研究.我正在评估一种头发陷阱方法与相机陷阱在检测3种(红松鼠,松貂和入侵灰松鼠)方面的成功.我想看看是什么影响了各种物种的检测(或非检测).一个假设是在现场检测到另一个焦点物种会影响红松鼠的可探测性.鉴于松貂是红松鼠的捕食者并且灰松鼠是竞争者,这两个物种在一个地点的存在可能会影响红松鼠的可探测性.

这会显示概率吗？ inv.logit(-1.14 - 0.1322 * nonRS events)

 glm(formula = RS_sticky ~ NonRSevents_before1stRS, family = binomial(link = "logit"), data = data)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7432  -0.7432  -0.7222  -0.3739   2.0361  
Coefficients:
                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)              -1.1455     0.4677  -2.449   0.0143 *
NonRSevents_before1stRS  -0.1322     0.1658  -0.797   0.4255  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
   (Dispersion parameter for binomial family taken to be 1)
Null deviance: 34.575  on 33  degrees of freedom
   Residual deviance: 33.736  on 32  degrees of freedom
  (1 observation deleted due to missingness)
   AIC: 37.736
  Number of Fisher Scoring iterations: 5*

Answer 1

如果要预测预测变量的指定值集的响应概率:

pframe <- data.frame(NonRSevents_before1stRS=4)
predict(fitted_model, newdata=pframe, type="response")

fitted_model你的glm()拟合结果在哪里,你存储在一个变量中.你可能不熟悉使用R的方法来进行统计分析,这对拟合模型存储为对象/在一个变量,然后应用不同的方法来它(summary(),plot(),predict(),residuals(),...)

• 这显然只是一个简单的例子:我不知道4是否是NonRSevents_before1stRS变量的合理值)
• 您可以指定更多不同的值来同时进行预测(data.frame(NonRSevents_before1stRS=c(4,5,6,7,8)))
• 如果您有多个预测变量,则必须为每个预测指定一些值,例如 data.frame(x=4:8,y=mean(orig_data$y), ...)

如果您想要原始数据集中观测值的预测概率,那么就是 predict(fitted_model, type="response")

你是正确的inv.logit()(从一堆不同的包,不知道你正在使用哪个)或plogis()(从基础R,基本上相同)将从logit或log-odds量表转换为概率量表,所以

plogis(predict(fitted_model))

也可以工作(predict默认情况下提供关于链接函数的预测[在这种情况下是logit/log-odds]比例).

Answer 2

逻辑回归中的因变量是对数优势比。我们将说明如何使用MASS包裹中的航天飞机自动着陆器数据来解释系数。

加载数据后，我们将创建一个二进制因变量，其中：

1 = autolander used, 
0 = autolander not used. 

我们还将为航天飞机稳定性创建一个二元自变量：

1 = stable positioning
0 = unstable positioning. 

然后，我们将运行glm()与family=binomial(link="logit")。由于系数是对数优势比，我们将对它们取指数以将它们变回优势比。

library(MASS)
str(shuttle)
shuttle$stable <- 0
shuttle[shuttle$stability =="stab","stable"] <- 1
shuttle$auto <- 0
shuttle[shuttle$use =="auto","auto"] <- 1

fit <- glm(use ~ factor(stable),family=binomial(link = "logit"),data=shuttle) # specifies base as unstable

summary(fit)
exp(fit$coefficients)

...和输出：

> fit <- glm(use ~ factor(stable),family=binomial(link = "logit"),data=shuttle) # specifies base as unstable
> 
> summary(fit)

Call:
glm(formula = use ~ factor(stable), family = binomial(link = "logit"), 
data = shuttle)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1774  -1.0118  -0.9566   1.1774   1.4155  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)  
(Intercept)      4.747e-15  1.768e-01   0.000   1.0000  
factor(stable)1 -5.443e-01  2.547e-01  -2.137   0.0326 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 350.36  on 255  degrees of freedom
Residual deviance: 345.75  on 254  degrees of freedom
AIC: 349.75

Number of Fisher Scoring iterations: 4

> exp(fit$coefficients)
    (Intercept) factor(stable)1 
      1.0000000       0.5802469 
> 

0 的截距是不稳定的对数几率，系数 -.5443 是稳定的对数几率。在对系数求幂后，我们观察到在航天飞机不稳定的情况下使用自动着陆器的几率为 1.0，如果航天飞机稳定，则乘以 0.58。这意味着如果航天飞机具有稳定的定位，则不太可能使用自动着陆器。

计算自动着陆器使用的概率

我们可以通过两种方式做到这一点。首先，手动方法：使用以下等式对系数求幂并将几率转换为概率。

p = odds / (1 + odds) 

使用航天飞机自动着陆器数据，它的工作方式如下。

# convert intercept to probability
odds_i <- exp(fit$coefficients[1])
odds_i / (1 + odds_i)
# convert stable="stable" to probability
odds_p <- exp(fit$coefficients[1]) * exp(fit$coefficients[2])
odds_p / (1 + odds_p)

...和输出：

> # convert intercept to probability
> odds_i <- exp(fit$coefficients[1])
> odds_i / (1 + odds_i)
(Intercept) 
        0.5 
> # convert stable="stable" to probability
> odds_p <- exp(fit$coefficients[1]) * exp(fit$coefficients[2])
> odds_p / (1 + odds_p)
(Intercept) 
  0.3671875 
>

航天飞机不稳定时使用自动着陆器的概率为 0.5，航天飞机稳定时使用自动着陆器的概率为 0.37。

生成概率的第二种方法是使用predict()函数。

# convert to probabilities with the predict() function
predict(fit,data.frame(stable="0"),type="response")
predict(fit,data.frame(stable="1"),type="response")

请注意，输出与手动计算的概率匹配。

> # convert to probabilities with the predict() function
> predict(fit,data.frame(stable="0"),type="response")
  1 
0.5 
> predict(fit,data.frame(stable="1"),type="response")
        1 
0.3671875 
> 

将此应用于 OP 数据

我们可以将这些步骤应用于glm()OP的输出，如下所示。

coefficients <- c(-1.1455,-0.1322)
exp(coefficients)
odds_i <- exp(coefficients[1])
odds_i / (1 + odds_i)
# convert nonRSEvents = 1 to probability
odds_p <- exp(coefficients[1]) * exp(coefficients[2])
odds_p / (1 + odds_p)
# simulate up to 10 nonRSEvents prior to RS
coef_df <- data.frame(nonRSEvents=0:10,
                  intercept=rep(-1.1455,11),
                  nonRSEventSlope=rep(-0.1322,11))
coef_df$nonRSEventValue <- coef_df$nonRSEventSlope * 
coef_df$nonRSEvents
coef_df$intercept_exp <- exp(coef_df$intercept)
coef_df$slope_exp <- exp(coef_df$nonRSEventValue)
coef_df$odds <- coef_df$intercept_exp * coef_df$slope_exp
coef_df$probability <- coef_df$odds / (1 + coef_df$odds)
# print the odds & probabilities by number of nonRSEvents
coef_df[,c(1,7:8)]

...以及最终输出。

> coef_df[,c(1,7:8)]
   nonRSEvents    odds probability
1            0 0.31806     0.24131
2            1 0.27868     0.21794
3            2 0.24417     0.19625
4            3 0.21393     0.17623
5            4 0.18744     0.15785
6            5 0.16423     0.14106
7            6 0.14389     0.12579
8            7 0.12607     0.11196
9            8 0.11046     0.09947
10           9 0.09678     0.08824
11          10 0.08480     0.07817
>