| Distribuicao tolerancia | Ligacao | Forma |
|---|---|---|
| Uniforme | Identidade | Linear restrita |
| Normal | Probit | Simetrica (caudas leves) |
| Logistica | Logit | Simetrica (caudas pesadas) |
| Valor extremo | Log-log compl. | Assimetrica |
Modelos Lineares Generalizados
Apostila — Capítulo 7: Variáveis Binárias e Regressão Logística
Nota ao leitor. Esta apostila resume o Capítulo 7 de An Introduction to Generalized Linear Models (Dobson & Barnett, 4ª ed., 2018). O capítulo cobre Regressão Logística e modelos dose-resposta para respostas binárias e binomiais. Inclui modelos probit, logístico e log-log complementar; inferência, resíduos, diagnósticos e interpretação de razões de chances vs. razões de prevalência.
Este material foi produzido em Quarto Markdown, com apoio de ferramentas de inteligência artificial na organização e síntese do conteúdo, tendo sido integralmente revisado e validado pelo autor.
1 Distribuições de Probabilidade para Respostas Binárias
1.1 Variável Bernoulli
Em muitos estudos o desfecho assume apenas dois valores — vivo / morto, presente / ausente. Chamamos genericamente de sucesso e falha. A peça fundamental é a variável de Bernoulli:
\[ Z = \begin{cases} 1 & \text{sucesso} \\ 0 & \text{fracasso} \end{cases} \]
com \(\Pr(Z=1)=\pi\). Distribuição de Bernoulli \(B(\pi)\).
1.2 Distribuição Binomial
Somando n ensaios independentes com a mesma probabilidade ππ, o número de sucessos \(Y = \sum Z_{j}\) segue modelo binomial Bin(n,\(\pi\))
\[ \Pr(Y=y)=\binom{n}{y}\pi^{y}(1-\pi)^{n-y}, \quad y=0,1,\ldots,n \tag{7.2} \]
\(E(Y)=n\pi\), \(\text{var}(Y)=n\pi(1-\pi)\).
por que isso importa
A verossimilhança conjunta de N subgrupos binomiais é membro da família exponencial — é exatamente o que permite tratar dados binários como um GLM e reaproveitar toda a maquinaria dos Capítulos 3–5.
1.3 \(N\) Grupos Independentes
Se \(Y_{i}\sim\text{Bin}(n_{i},\pi_{i})\), a log-verossimilhança é:
\[ \ell(\boldsymbol{\pi};\mathbf{y})=\sum_{i=1}^{N}\left[y_{i}\log\frac{\pi_{i}}{1-\pi_{i}}+n_{i}\log(1-\pi_{i})+\log\binom{n_{i}}{y_{i}}\right] \tag{7.3} \]
2 MLGs para Respostas Binárias
2.1 Função de Ligação
\[g(\pi_{i})=\mathbf{x}_{i}^{T}\boldsymbol{\beta}\]
logit
\[\log \dfrac{\pi}{1-\pi} = \beta_{1} + \beta_{2}x\] Distribuição logística. O termo \(\log \dfrac{\pi}{1-\pi}\) é o log das chances.
probit
\[\Phi^{-1}(\pi) = \beta_{1}+\beta_{2}x\] Distribuição Normal. O termo \(x=\mu\) é a dose letal mediana LD(50).
cloglog
\[\log \{-\log(1-\pi)\} = \beta_{1}+\beta_{2}x\] Valor extremo. Assimétrico — difere nas caudas (\(\pi\) perto de 0 ou 1).
Código
x_seq <- seq(-4, 4, length.out = 300)
df_lk <- data.frame(
x = rep(x_seq, 3),
pi = c(plogis(x_seq), pnorm(x_seq), 1-exp(-exp(x_seq))),
modelo = rep(c("Logistico","Probit","Log-log compl."), each=300)
)
ggplot(df_lk, aes(x=x, y=pi, color=modelo, linetype=modelo)) +
geom_line(linewidth=1.1) +
scale_color_manual(values=c("#4C72B0","#C44E52","#55A868"), name="Modelo") +
scale_linetype_manual(values=c("solid","dashed","dotdash"), name="Modelo") +
labs(x="Preditor linear", y=expression(pi(x)),
title="Curvas dose-resposta") + tema
3 Modelos Dose-Resposta
3.1 Exemplo 7.3.1 — Mortalidade de Besouros
Bliss (1935) expôs besouros a dissulfeto de carbono gasoso em 8 concentrações. Para cada dose \(x_{i}\) (em \(\log_{10}\)) temos \(n_{i}\) besouros e \(y_{i}\) mortos.
Código
x_b <- c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)
n_b <- c(59,60,62,56,63,59,62,60)
y_b <- c(6,13,18,28,52,53,61,60)
p_b <- y_b/n_b
dados_b <- data.frame(x=x_b, n=n_b, y=y_b, p=p_b)| Dose xi | Total ni | Mortos yi | Proporcao pi |
|---|---|---|---|
| 1.6907 | 59 | 6 | 0.1017 |
| 1.7242 | 60 | 13 | 0.2167 |
| 1.7552 | 62 | 18 | 0.2903 |
| 1.7842 | 56 | 28 | 0.5000 |
| 1.8113 | 63 | 52 | 0.8254 |
| 1.8369 | 59 | 53 | 0.8983 |
| 1.8610 | 62 | 61 | 0.9839 |
| 1.8839 | 60 | 60 | 1.0000 |
Código
ggplot(dados_b, aes(x=x, y=p)) +
geom_point(size=3.5, color="blue4", alpha=0.9) +
scale_y_continuous(limits=c(0,1), labels=scales::percent) +
labs(x="Dose xi", y="Proporcao morta",
title="Mortalidade de besouros vs. dose") + tema
Código
mat_b <- cbind(y_b, n_b-y_b)
mod_logit <- glm(mat_b~x_b, family=binomial(link="logit"))
mod_probit <- glm(mat_b~x_b, family=binomial(link="probit"))
mod_cloglog<- glm(mat_b~x_b, family=binomial(link="cloglog"))
yhat_logit <- fitted(mod_logit) * n_b
yhat_probit <- fitted(mod_probit) * n_b
yhat_cloglog <- fitted(mod_cloglog) * n_b
D_logit <- deviance(mod_logit)
D_probit <- deviance(mod_probit)
D_cloglog <- deviance(mod_cloglog)
cf_l <- coef(summary(mod_logit))
cf_p <- coef(summary(mod_probit))
cf_c <- coef(summary(mod_cloglog))| Observado | Logistico | Probit | Log-log compl. |
|---|---|---|---|
| 6 | 3.46 | 3.36 | 5.59 |
| 13 | 9.84 | 10.72 | 11.28 |
| 18 | 22.45 | 23.48 | 20.95 |
| 28 | 33.90 | 33.82 | 30.37 |
| 52 | 50.10 | 49.62 | 47.78 |
| 53 | 53.29 | 53.32 | 54.14 |
| 61 | 59.22 | 59.66 | 61.11 |
| 60 | 58.74 | 59.23 | 59.95 |
| Parametro | Logistico | Probit | Log-log compl. |
|---|---|---|---|
| b1 (intercepto) | -60.72 (5.18) | -34.94 (2.65) | -39.57 (3.24) |
| b2 (dose) | 34.27 (2.91) | 19.73 (1.49) | 22.04 (1.80) |
| Deviance D | 11.23 | 10.12 | 3.45 |
| gl (N-p) | 6 | 6 | 6 |
Código
x_seq2 <- seq(min(x_b)-0.02, max(x_b)+0.02, length.out=300)
df_fit <- data.frame(
x = rep(x_seq2,3),
pi = c(predict(mod_logit, newdata=data.frame(x_b=x_seq2), type="response"),
predict(mod_probit, newdata=data.frame(x_b=x_seq2), type="response"),
predict(mod_cloglog,newdata=data.frame(x_b=x_seq2), type="response")),
modelo=rep(c(sprintf("Logistico (D=%.2f)",D_logit),
sprintf("Probit (D=%.2f)",D_probit),
sprintf("Log-log compl. (D=%.2f)",D_cloglog)), each=300)
)
ggplot() +
geom_point(data=dados_b, aes(x=x,y=p), size=3.5, color="gray30", alpha=0.8) +
geom_line(data=df_fit, aes(x=x,y=pi,color=modelo,linetype=modelo), linewidth=1.05) +
scale_color_manual(values=c("blue4","red4","green4"), name="Modelo") +
scale_linetype_manual(values=c("solid","dashed","dotdash"), name="Modelo") +
scale_y_continuous(limits=c(0,1), labels=scales::percent) +
labs(x="Dose xi", y="Proporcao morta", title="Besouros: curvas ajustadas") + tema
Conclusao: Log-log complementar ajusta melhor (\(D=3{,}45\) vs. \(D\approx11\)).
4 Modelo Logistico Geral
4.1 Definicao
A reta simples \(\beta_{1}+\beta_{2}x\) é só um caso particular. No modelo geral, o logit é uma combinação linear de covariáveis e variáveis indicadoras de fatores:
\[ \text{logit}(\pi_{i})=\log\frac{\pi_{i}}{1-\pi_{i}}=\mathbf{x}_{i}^{T}\boldsymbol{\beta} \tag{7.4a} \]
\[ \pi_{i}=\frac{e^{\mathbf{x}_{i}^{T}\boldsymbol{\beta}}}{1+e^{\mathbf{x}_{i}^{T}\boldsymbol{\beta}}} \tag{7.4b} \]
4.2 Deviance Binomial
O modelo geral éo análogo da regressão múltipla e da ANOVA, agora para respostas binárias. A estimação é a mesma quer os dados venham agrupados (contagens \(y_{i}\) em \(n_{i}\) ensaios) quer binários (cada linha 0/1, com \(n_{i}=1\)). A deviance não envolve parâmetros de perturbação, então testamos hipóteses diretamente com
\[ D=2\sum_{i=1}^{N}\left[y_{i}\log\frac{y_{i}}{\hat{y}_{i}}+(n_{i}-y_{i})\log\frac{n_{i}-y_{i}}{n_{i}-\hat{y}_{i}}\right]\sim\chi^{2}(N-p) \tag{7.5} \]
5 Exemplo 7.4.1 — Anteras Embriogenicas
Anteras de Datura innoxia sob tratamento (armazenamento a 3\(^{o}\)C) vs. controle, em 3 forças de centrifugação.
Código
y_a <- c(55,52,57,55,50,50)
n_a <- c(102,99,108,76,81,90)
x_k <- c(log(40),log(150),log(350))
x_a <- rep(x_k,2)
stor <- factor(rep(c("Controle","Tratamento"), each=3))
dados_a <- data.frame(y=y_a, n=n_a, x=x_a, storage=stor, p=y_a/n_a)| Condicao | xk=ln(g) | yjk | njk | pjk |
|---|---|---|---|---|
| Controle | 3.689 | 55 | 102 | 0.539 |
| Controle | 5.011 | 52 | 99 | 0.525 |
| Controle | 5.858 | 57 | 108 | 0.528 |
| Tratamento | 3.689 | 55 | 76 | 0.724 |
| Tratamento | 5.011 | 50 | 81 | 0.617 |
| Tratamento | 5.858 | 50 | 90 | 0.556 |
Código
ggplot(dados_a, aes(x=x, y=p, color=storage, shape=storage)) +
geom_point(size=4, alpha=0.9) +
scale_color_manual(values=c("#4C72B0","#C44E52"), name="Condicao") +
scale_shape_manual(values=c(16,17), name="Condicao") +
scale_y_continuous(limits=c(0.45,0.80), labels=scales::percent) +
labs(x="xk = ln(g)", y="Proporcao embriogenica",
title="Anteras embriogenicas") + tema
Código
mat_a <- cbind(y_a, n_a-y_a)
newstor <- as.numeric(stor)-1
mod1 <- glm(mat_a~newstor+x_a+newstor:x_a, family=binomial(link="logit"))
mod2 <- glm(mat_a~newstor+x_a, family=binomial(link="logit"))
mod3 <- glm(mat_a~x_a, family=binomial(link="logit"))
cf1 <- coef(summary(mod1))
cf2 <- coef(summary(mod2))
cf3 <- coef(summary(mod3))Comparamos três modelos para o logit, do mais flexível ao mais simples. A diferença de deviances entre modelos aninhados é o teste de razão de verossimilhança:
| Modelo | a1 (EP) | a2-a1 (EP) | b (EP) | b2-b1 (EP) | D | gl |
|---|---|---|---|---|---|---|
| Modelo 1: aj + bj*xk | 0.234 (0.628) | 1.977 (0.998) | -0.023 (0.127) | -0.319 (0.199) | 0.028 | 2 |
| Modelo 2: aj + b*xk | 0.877 (0.487) | 0.407 (0.175) | -0.155 (0.097) | --- | 2.619 | 3 |
| Modelo 3: a + b*xk | 1.021 (0.481) | --- | -0.148 (0.096) | --- | 8.092 | 4 |
Código
pv_slopes <- pchisq(deviance(mod2)-deviance(mod1), df=1, lower.tail=FALSE)
pv_stor <- pchisq(deviance(mod3)-deviance(mod2), df=1, lower.tail=FALSE)Inclinacoes iguais (Mod 1 vs 2): \(\Delta D = 2.591\), \(p = 0.107\) — sem evidencia contra \(H_0\).
Efeito armazenamento (Mod 2 vs 3): \(\Delta D = 5.473\), \(p = 0.019\) — efeito significativo.
| Condicao | xk | Observado | Modelo 1 | Modelo 2 | Modelo 3 |
|---|---|---|---|---|---|
| Controle | 3.69 | 55 | 54.82 | 58.75 | 62.91 |
| Controle | 5.01 | 52 | 52.47 | 52.03 | 56.40 |
| Controle | 5.86 | 57 | 56.72 | 53.22 | 58.18 |
| Tratamento | 3.69 | 55 | 54.83 | 51.01 | 46.88 |
| Tratamento | 5.01 | 50 | 50.43 | 50.59 | 46.14 |
| Tratamento | 5.86 | 50 | 49.74 | 53.40 | 48.49 |
6 Estatisticas de Bondade de Ajuste
6.1 Pearson e Deviance
\[X^{2}=\sum_{k=1}^{m}\frac{(y_{k}-n_{k}\hat{\pi}_{k})^{2}}{n_{k}\hat{\pi}_{k}(1-\hat{\pi}_{k})} \tag{7.6}\]
Ambas \(D\) e \(X^{2}\) \(\sim\chi^{2}(N-p)\).
6.1.1 Por que \(D\approx X^{2}\)?
Expansao de Taylor de \(s\log(s/t)\) em torno de \(s=t\):
\[s\log(s/t)=(s-t)+\frac{(s-t)^{2}}{2t}+\ldots\]
Os termos lineares se cancelam pelas equacoes de score; os quadraticos dao:
\[D\approx\sum_{k}\frac{(y_{k}-n_{k}\hat{\pi}_{k})^{2}}{n_{k}\hat{\pi}_{k}(1-\hat{\pi}_{k})}=X^{2}\]
6.2 AIC e BIC
\[\text{AIC}=-2\,\ell(\hat{\boldsymbol{\pi}};\mathbf{y})+2p \tag{7.7}\]
| Estatistica | Formula | Distrib. sob H0 | Uso |
|---|---|---|---|
| D (Deviance) | 2 sum oi*log(oi/ei) | chi2(N-p) | Bondade ajuste |
| X2 (Pearson) | sum(oi-ei)^2/ei | chi2(N-p) | Bondade ajuste |
| C (vs. minimo) | D0 - D1 | chi2(p-1) | Testa preditores |
| Pseudo-R2 | 1 - l(pi_hat)/l(pi_tilde) | --- | Melhoria log-lik |
| AIC | -2l + 2p | --- | Comparacao modelos |
| BIC | -2l + p*ln(N) | --- | Comparacao (penaliza) |
7 Residuos e Diagnosticos
7.1 Tipos de Residuos
Pearson: \(X_{k}=(y_{k}-n_{k}\hat{\pi}_{k})/\sqrt{n_{k}\hat{\pi}_{k}(1-\hat{\pi}_{k})}\) \(\quad(7.8)\)
Deviance: \(d_{k}=\text{sign}(y_{k}-n_{k}\hat{\pi}_{k})\sqrt{2[\ldots]}\) \(\quad(7.9)\)
Padronizados: \(r_{Pk}=X_{k}/\sqrt{1-h_{k}}\), \(\quad r_{Dk}=d_{k}/\sqrt{1-h_{k}}\)
7.2 Hosmer-Lemeshow
e quando os dados são binários?
Se cada observação tem covariáveis próprias (\(y_{i} \in \{0,1\}\)), nem D nem \(X^{2}\) servem. A saída de Hosmer-Lemeshow : agrupar as observações em \(\sim 10\) faixas pela probabilidade prevista e calcular um \(X^{2}\) de contingência \(g \times 2\), com \(X^{2}_{HL} \sim \chi^{2}_{g-2}\).
\[X^{2}_{HL}=\sum_{\text{celulas}}\frac{(o-e)^{2}}{e}\sim\chi^{2}(g-2)\]
8 Exemplo 7.8 — Senilidade e WAIS
- 54 idosos: sintomas de senilidade (sim/não) vs. escore WAIS.
Senilidade
Envelhecimento patológico (decorrente de doenças).
WAIS
Escala de Inteligência Wechsler para Adultos (do inglês Wechsler Adult Intelligence Scale ou WAIS) é uma medida padronizada (QI) e de índices cognitivos.
Código
x_s <- c(
9, 7, 7,17,13,
13, 5,16,14,13,
6,14, 9,19, 9,
8,13, 9, 9,15,
10,16,11,11,10,
4,10,13,14,11,
14,12,15,10,12,
8,11,13,16, 4,
11,14,10,10,14,
7,15,11,16,20,
9,18, 6,14
)
s_wais <- c(
1,1,0,0,0,
1,1,0,0,0,
1,1,0,0,0,
1,0,0,0,0,
1,0,0,0,0,
1,0,0,0,0,
1,0,0,0,0,
1,0,0,0,0,
1,0,0,0,0,
1,0,0,0,0,
1,0,0,0
)
stopifnot(length(x_s)==54, length(s_wais)==54, sum(s_wais)==14)
dados_s <- data.frame(x=x_s, s=s_wais)
library(dplyr)
wais_grp <- dados_s |>
group_by(x) |>
summarise(y=sum(s), n=n(), .groups="drop") |>
arrange(x)
wais_grp$p_obs <- wais_grp$y / wais_grp$nCódigo
mod_s <- glm(cbind(y,n-y)~x, family=binomial(link="logit"), data=wais_grp)
wais_grp$pi_hat <- fitted(mod_s)
cf_s <- coef(summary(mod_s))
wais_grp$X_k <- (wais_grp$y - wais_grp$n*wais_grp$pi_hat) /
sqrt(wais_grp$n*wais_grp$pi_hat*(1-wais_grp$pi_hat))
wais_grp$d_k <- residuals(mod_s, type="deviance")| x (WAIS) | yi | ni | pi_hat | Xk (Pearson) | dk (deviance) |
|---|---|---|---|---|---|
| 4 | 1 | 2 | 0.752 | -0.826 | -0.766 |
| 5 | 1 | 1 | 0.687 | 0.675 | 0.866 |
| 6 | 1 | 2 | 0.614 | -0.330 | -0.326 |
| 7 | 2 | 3 | 0.535 | 0.458 | 0.464 |
| 8 | 2 | 2 | 0.454 | 1.551 | 1.777 |
| 9 | 2 | 6 | 0.376 | -0.214 | -0.216 |
| 10 | 1 | 6 | 0.303 | -0.728 | -0.771 |
| 11 | 1 | 6 | 0.240 | -0.419 | -0.436 |
| 12 | 0 | 2 | 0.186 | -0.675 | -0.906 |
| 13 | 1 | 6 | 0.142 | 0.176 | 0.172 |
| 14 | 2 | 7 | 0.107 | 1.535 | 1.306 |
| 15 | 0 | 3 | 0.080 | -0.509 | -0.705 |
| 16 | 0 | 4 | 0.059 | -0.500 | -0.696 |
| 17 | 0 | 1 | 0.043 | -0.213 | -0.297 |
| 18 | 0 | 1 | 0.032 | -0.181 | -0.254 |
| 19 | 0 | 1 | 0.023 | -0.154 | -0.216 |
| 20 | 0 | 1 | 0.017 | -0.131 | -0.184 |
| Soma dos quadrados: X2 = 8.083, D = 9.419 |
Código
x_seq3 <- seq(min(x_s)-0.5, max(x_s)+0.5, length.out=300)
pi_seq <- predict(mod_s, newdata=data.frame(x=x_seq3), type="response")
ggplot() +
geom_point(data=wais_grp, aes(x=x, y=p_obs, size=n),
color="blue4", alpha=0.75) +
geom_line(data=data.frame(x=x_seq3,pi=pi_seq),
aes(x=x,y=pi), color="red4", linewidth=1.1, linetype="dashed") +
scale_size_continuous(range=c(2,8), name="n") +
scale_y_continuous(limits=c(-0.05,1.05), labels=scales::percent) +
labs(x="Pontuacao WAIS", y="Proporcao com sintomas",
title="Senilidade e WAIS",
subtitle=paste0("b1=",round(cf_s[1,1],3),
" (EP=",round(cf_s[1,2],3),") ",
"b2=",round(cf_s[2,1],4),
" (EP=",round(cf_s[2,2],4),")")) + tema
Código
X2_w <- sum(wais_grp$X_k^2)
D_w <- deviance(mod_s)
mod_s0 <- glm(cbind(y,n-y)~1, family=binomial(link="logit"), data=wais_grp)
C_w <- as.numeric(2*(logLik(mod_s)-logLik(mod_s0)))
pC_w <- pchisq(C_w, df=1, lower.tail=FALSE)
wais_grp$grupo_hl <- cut(wais_grp$pi_hat,
breaks=c(0,0.107,0.303,1), labels=c("ate 0.107","0.108 a 0.303","acima 0.303"))
hl_tab <- wais_grp |>
group_by(grupo_hl) |>
summarise(o_sim=sum(y), e_sim=round(sum(n*pi_hat),3),
o_nao=sum(n-y), e_nao=round(sum(n*(1-pi_hat)),3),
n_tot=sum(n), .groups="drop")
X2_HL <- sum((hl_tab$o_sim-hl_tab$e_sim)^2/hl_tab$e_sim +
(hl_tab$o_nao-hl_tab$e_nao)^2/hl_tab$e_nao)
pHL <- pchisq(X2_HL, df=nrow(hl_tab)-2, lower.tail=FALSE)
pseudo_r2 <- as.numeric(1-logLik(mod_s)/logLik(mod_s0))
aic_w <- AIC(mod_s)| Faixa pi_hat | Com sint.obs | Com sint.esp | Sem sint.obs | Sem sint.esp | Total |
|---|---|---|---|---|---|
| ate 0.107 | 2 | 1.335 | 16 | 16.665 | 18 |
| 0.108 a 0.303 | 2 | 2.659 | 12 | 11.341 | 14 |
| acima 0.303 | 10 | 10.006 | 12 | 11.994 | 22 |
| X2_HL = 0.56 ~ chi2(1); p-valor = 0.454 |
| Estatistica | Valor |
|---|---|
| b1 (intercepto) | 2.404 (EP=1.192) |
| b2 (inclinacao) | -0.3235 (EP=0.1140) |
| X2 Pearson (gl=15) | 8.083 (p=0.920) |
| D deviance (gl=15) | 9.419 (p=0.855) |
| C vs. minimo (gl=1) | 10.789 (p<0.001) |
| Pseudo-R2 | 0.3120 |
| AIC | 27.792 |
| X2_HL Hosmer-Lemeshow (gl=1) | 0.56 (p=0.454) |
9 Razoes de Chances e Razoes de Prevalencia
Com a ligação logit, os coeficientes exponenciados são razões de chances (odds ratios) — e são constantemente mal interpretados como riscos relativos. A razão de chances compara as chances \(\pi/(1-\pi)\); a razão de prevalência (risco relativo) compara as probabilidades.
Para dois grupos com probabilidades \(\pi_{1}\) e \(\pi_{2}\):
\[\text{RC}=\frac{\pi_{1}/(1-\pi_{1})}{\pi_{2}/(1-\pi_{2})}, \qquad \text{RP}=\frac{\pi_{1}}{\pi_{2}}\]
No modelo logistico \(g(\pi_{i})=\beta_{1}+\beta_{2}x_{i}\):
\[\text{RC}=e^{\beta_{2}}\]
- Modelo logístico ajustado
\[\log \dfrac{\pi}{1-\pi} = 2,404 - 0,324 \times (WAIS)\]
A variável explicativa é a pontuação no WAIS (inteiro de 4 a 20) e a resposta é ter sintomas de senilidade (\(\pi = 1\)) ou não.
OR: O coeficiente \(b_{2} = -0,324\) é o log-odds ratio associado a um aumento de 1 ponto no WAIS. Portanto:
\[OR = e^{b_{2}} = e^{-0,324} \approx 0,723\]
Interpretação
Cada ponto a mais na pontuação WAIS multiplica as chances de ter sintomas de senilidade por 0,723 — uma redução de 27,7%.”
Às vezes, 1 ponto no WAIS é uma mudança pequena. Podemos interpretar uma mudança de 5 pontos:
\[OR = e^{5\times b_{2}} = e^{5\times(-0,324)} \approx 0,198\]
O que isso significa
Pessoas com 5 pontos a mais no WAIS têm aproximadamente 0.198 vezes a chance de apresentar sintomas de senilidade. Isto é, um aumento de 5 pontos no WAIS está associado a uma redução de cerca de 80,2% (1-0,198) nas chances de sintomas.
- No sentido inverso: (\(1/0,198\approx 5,04\))
Pessoas com 5 pontos a menos no WAIS têm cerca de 5 vezes a chance de apresentar sintomas de senilidade.
A RP compara a probabilidade de sintomas de senilidade entre dois valores do WAIS. Diferentemente do OR, a RP não compara chances. Ela compara probabilidades. Por exempo, usando o modelo ajustado:
\[ \hat{\pi}(x) = \dfrac{e^{2,404-0,3235 x}}{1+e^{2,404-0,3235 x}} \] Comparando WAIS 5 com WAIS 10:
\[RP = \dfrac{\hat{\pi}(10)}{\hat{\pi}(5)} = \dfrac{0,303}{0,687} \approx 0,442\]
Interpretação
Uma pessoa com WAIS 10 tem aproximadamente 44% da probabilidade de sintomas de uma pessoa com WAIS 5. Ou, de forma equivalente, a probabilidade ajustada de sintomas é cerca de 55,8% menor em WAIS 10 do que em WAIS 5.
Resumindo
- No exemplo Senilidade e WAIS, o OR mostra como o escore WAIS altera as chances de sintomas. Como OR<1, maiores escores WAIS estão associados a menores chances de senilidade. Já a RP mostra como o WAIS altera diretamente a probabilidade de sintomas. A RP costuma ser mais intuitiva, mas no modelo logístico o resultado natural é o OR.
- OR é a razão entre chances; RP é a razão entre probabilidades. Quando o evento é comum, como sintomas de senilidade em baixos valores de WAIS, OR e RP podem ser bem diferentes.
Código
WAIS=c(5,10,15,10)
prob.ajust <- (exp(2.404-0.3235*WAIS))/(1+exp(2.404-0.3235*WAIS))
tab_RP = data.frame(WAIS,round(prob.ajust,3))
kable(tab_RP,
caption="**Algumas probabilidade ajustadas",
col.names=c("WAIS(x)","Probabilidade ajustada"),
align="c") |>
kable_styling(bootstrap_options=c("striped","hover","condensed"),
full_width=FALSE, font_size=12)| WAIS(x) | Probabilidade ajustada |
|---|---|
| 5 | 0.687 |
| 10 | 0.303 |
| 15 | 0.080 |
| 10 | 0.303 |
Código
y_sum <- c(164,155); n_sum <- c(309,247); stor2 <- c(0,1)
mat_sum <- cbind(y_sum, n_sum-y_sum)
mod_logit2 <- glm(mat_sum~stor2, family=binomial(link="logit"))
mod_log2 <- glm(mat_sum~stor2, family=binomial(link="log"))| Condicao | Sucesso | Fracasso | Total |
|---|---|---|---|
| Controle | 164 | 145 | 309 |
| Tratamento | 155 | 92 | 247 |
| Total | 319 | 237 | 556 |
Código
cf_logit2 <- coef(summary(mod_logit2))
cf_log2 <- coef(summary(mod_log2))
OR_ci <- exp(confint.default(mod_logit2)[2,])
PR_ci <- exp(confint.default(mod_log2)[2,])
tab_orpr <- data.frame(
Modelo =c("Logistico (logit)","Log-binomial (log)"),
Medida =c("Razao de chances (OR)","Razao de prevalencia (RP)"),
Estimativa=c(round(exp(cf_logit2[2,1]),3), round(exp(cf_log2[2,1]),3)),
IC95 =c(sprintf("(%.3f, %.3f)",OR_ci[1],OR_ci[2]),
sprintf("(%.3f, %.3f)",PR_ci[1],PR_ci[2]))
)| Modelo | Medida | Estimativa | IC 95% |
|---|---|---|---|
| Logistico (logit) | Razao de chances (OR) | 1.490 | (1.059, 2.095) |
| Log-binomial (log) | Razao de prevalencia (RP) | 1.182 | (1.026, 1.363) |
Código
pi2 <- seq(0.05,0.65,length.out=200)
pi1 <- 1.5*pi2; pi1 <- pi1[pi1<1]; pi2 <- pi2[length(pi1)]
pi2_v <- seq(0.05,0.65,length.out=200)
pi1_v <- 1.5*pi2_v
ok <- pi1_v < 1
pi2_v <- pi2_v[ok]; pi1_v <- pi1_v[ok]
OR_v <- (pi1_v/(1-pi1_v))/(pi2_v/(1-pi2_v))
ggplot(data.frame(pi2=pi2_v,OR=OR_v), aes(x=pi2,y=OR)) +
geom_line(color="lightblue", linewidth=1.2) +
geom_hline(yintercept=1.5, color="red", linetype="dashed") +
annotate("text",x=0.55,y=1.65,label="RP = 1.5",color="red",size=3.5) +
labs(x="Probabilidade no grupo referencia (pi2)",
y="Razao de chances (OR)",
title="OR vs. RP: divergem com pi2 elevado") + tema
10 Resumo do Capitulo 7
Respostas binárias → Bernoulli/Binomial, dentro da família exponencial.
Uma ligação (logit, probit, cloglog) comprime a reta para [0, 1].
Estimamos \(\beta\) por máxima verossimilhança, iterativamente.
Avaliamos o ajuste com \(D\), \(X^{2}\), Hosmer–Lemeshow e resíduos.
Interpretamos com cuidado: odds ratio \(\neq\) risco relativo.
| Elemento | Resultado/Formula |
|---|---|
| Modelo | logit(pi_i) = xi' beta |
| Distribuicao | Yi ~ Bin(ni, pi_i) |
| Ligacao canonica | log[pi/(1-pi)] (logit) |
| Deviance | D = 2 sum[yi*log(yi/y_hat_i) + ...] |
| Escore U1 | U1 = sum(yi - ni*pi_hat_i) |
| Peso IRLS wi | wi = ni*pi_hat_i*(1-pi_hat_i) |
| Resp. ajustada zi | zi = eta_hat_i + (yi - ni*pi_hat_i)/wi |
| Diagnosticos | Pearson Xk, deviance dk, padronizados, H-L |
| OR vs RP | OR = exp(bj); RP via link log |
10.1 Referencias
Dobson, A. J. & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4 ed.). CRC Press.
Bliss, C. I. (1935). The calculation of the dosage-mortality curve. Annals of Applied Biology, 22(1), 134–167.
Hosmer, D. W. & Lemeshow, S. (1980). A goodness-of-fit test for the multiple logistic regression model. Communications in Statistics, 10(10), 1043–1069.