Questão 1
a)
Com o apoio da análise posterior citada no enunciado, a suposição de independência entre os lançamentos possibilita a construção de um modelo mais simples e facilmente interpretável. Todavia, se os lançamentos, por algum motivo, forem dependentes, o ajuste logístico não refletirá a probabilidade de falha real dos O.rings. Se esse for o caso, seria interessante conduzir novos estudos, procurando por possíveis correlações entre as características de cada um deles, sempre com muita cautela.
b)
challenger <- read.csv(file = "http://leg.ufpr.br/~lucambio/CE073/20222S/challenger.csv")
ajuste1 <- glm(cbind(O.ring, Number-O.ring) ~ Temp + Pressure, family = binomial(link = "logit"), data = challenger)
coef(ajuste1)## (Intercept) Temp Pressure
## 2.520194641 -0.098296750 0.008484021Com isso,
\[2,520 - 0,098*Temp + 0,008*Pressure = log \left(\frac{\hat\pi}{1-\hat\pi}\right)\]
c)
Para verificar a significância das variáveis Temp e Pressure através do TRV, podemos utilizar a função Anova do pacote car.
library(car)
Anova(ajuste1, kind = "LR")## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(O.ring, Number - O.ring)
## LR Chisq Df Pr(>Chisq)
## Temp 5.1838 1 0.0228 *
## Pressure 1.5407 1 0.2145
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1d)
ajuste2 <- glm(cbind(O.ring, Number-O.ring) ~ Temp, family = binomial(link = "logit"), data = challenger)
AIC(ajuste1, ajuste2)## df AIC
## ajuste1 3 36.10589
## ajuste2 2 35.64654Com a remoção da variável Pressure, não significativa no primeiro modelo, temos um ajuste melhor (menor valor de AIC, critério que penaliza a inclusão de variáveis desnecessárias). Em contrapartida, temos um modelo com apenas uma variável explicativa, desconsiderando informações que poderiam trazer valor ao problema.
Questão 2
a)
round(coef(ajuste2), 3)## (Intercept) Temp
## 5.085 -0.116Com isso,
\[5,085 - 0.116*Temp = log \left(\frac{\hat\pi}{1 - \hat\pi}\right)\]
b)
# Gráfico 1
pred <- predict(ajuste2, newdata = data.frame(Temp = 31:81), type = "response")
plot(pred ~ c(31:81), data = challenger, col = "#ff5454",
xlim = c(31,81), ylim = c(0,1), xlab = "Temperatura",
ylab = "pi.hat", type = "l", lwd = 1.5)
# Gráfico 2 (como X~Bin(n,pi), E(X) = n*pi...)
plot(pred*6 ~ c(31:81), data = challenger, col = "steelblue",
xlim = c(31,81), xlab = "Temperatura",
ylab = "Número Esperado de Falhas", type = "l", lwd = 1.5)
c)
pred_link <- predict(ajuste2, newdata = data.frame(Temp = 31:81), se.fit = T)
df_pred <- data.frame(Temp = c(31:81),
pi.hat = pred,
LowerLimit = exp(pred_link$fit - 1.96*pred_link$se.fit) / (1 + exp(pred_link$fit - 1.96*pred_link$se.fit)),
UpperLimit = exp(pred_link$fit + 1.96*pred_link$se.fit) / (1 + exp(pred_link$fit + 1.96*pred_link$se.fit)))
plot(pred ~ c(31:81), data = challenger, col = "#ff5454",
xlim = c(31,81), ylim = c(0,1), xlab = "Temperatura",
ylab = "pi.hat", type = "l", lwd = 1.8)
lines(df_pred$Temp, df_pred$LowerLimit, col = "steelblue", lty = 2, lwd = 1.8)
lines(df_pred$Temp, df_pred$UpperLimit, col = "steelblue", lty = 2, lwd = 1.8) 
Acredito que as bandas de confiança sejam mais largas para temperaturas menores por conta de que não temos informações registradas sobre temperaturas mais baixas (menor valor do dataset = 58), gerando erros padrões mais altos.
d)
Utilizando os resultados da questão anterior:
head(df_pred, n = 1)## Temp pi.hat LowerLimit UpperLimit
## 1 31 0.8177744 0.1595947 0.9906587Suposições:
Linearidade entre a temperatura e o logito de \(\pi\).
Variância constante (dos erros)
e)
fit_model <- function(data, indices, temperatura){
sample_data <- data[indices, ]
model <- glm(cbind(O.ring, Number-O.ring) ~ Temp, family = binomial(link = "logit"), data = sample_data)
pred <- predict(model, newdata = data.frame(Temp = temperatura), type = "response")
return(pred)
}
library(boot)
# Temp = 31
temp31 <- boot(data = challenger, statistic = fit_model, R = 1000, temperatura = 31)
boot.ci(temp31, type = c("perc"), conf = 0.90)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = temp31, conf = 0.9, type = c("perc"))
##
## Intervals :
## Level Percentile
## 90% ( 0.1497, 0.9918 )
## Calculations and Intervals on Original Scale# Temp = 72
temp72 <- boot(data = challenger, statistic = fit_model, R = 1000, temperatura = 72)
boot.ci(temp72, type = c("perc"), conf = 0.90)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = temp72, conf = 0.9, type = c("perc"))
##
## Intervals :
## Level Percentile
## 90% ( 0.0057, 0.0780 )
## Calculations and Intervals on Original Scalef)
ajuste3 <- glm(cbind(O.ring, Number-O.ring) ~ Temp + I(Temp^2), family = binomial(link = "logit"), data = challenger)
AIC(ajuste2, ajuste3)## df AIC
## ajuste2 2 35.64654
## ajuste3 3 37.15184A inclusão de um termo quadrático resultou em um valor de AIC mais alto, implicando que a inclusão do mesmo não é necessária.
Questão 3
a)
ajuste_probit <- glm(cbind(O.ring, Number-O.ring) ~ Temp, family = binomial(link = "probit"), data = challenger)
ajuste_cloglog <- glm(cbind(O.ring, Number-O.ring) ~ Temp, family = binomial(link = "cloglog"), data = challenger)
pred_challenger <- data.frame(Temp = c(31:81),
PredLogit = predict(ajuste2, type = "response", newdata = data.frame(Temp = c(31:81))),
PredProbit = predict(ajuste_probit, type = "response", newdata = data.frame(Temp = c(31:81))),
PredClogLog = predict(ajuste_cloglog, type = "response", newdata = data.frame(Temp = c(31:81))))
# Gráfico
plot(pred_challenger$PredLogit ~ pred_challenger$Temp, type = "l", col = "steelblue", lwd = 2, ylim = c(0,1), xlab = "Temperatura", ylab = "pi.hat")
lines(pred_challenger$PredProbit ~ pred_challenger$Temp, col = "#ff5454", lwd = 2, ylim = c(0,1))
lines(pred_challenger$PredClogLog ~ pred_challenger$Temp, col = "orange", lwd = 2, ylim = c(0,1))
legend("topright", legend = c("Logit", "Probit", "Clog-log"), col = c("steelblue", "#ff5454", "orange"), lwd = 2)
Os ajustes probito e clog-log atuaram de forma muito similar no cálculo das probabilidades de falha para temperaturas mais altas, mas apresentaram um certo desvio para valores mais baixos.
b)
Decidi utilizar a função de ligação paramétrica Pregibin para a comparação.
library(glmx); library(gld)
pregibin <- function(shape) binomial(link = pregibon(shape[1], shape[2]))
ajuste_pregibin <- glmx(formula = cbind(O.ring, Number-O.ring) ~ Temp, data = challenger,
family = pregibin, xstart = c(0,0), xlink = "logit")
pred_challenger$PredPregibin <- predict(ajuste_pregibin, type = "response", newdata = data.frame(Temp = c(31:81)))
# Gráfico
par(mfrow = c(1,2))
plot(pred_challenger$PredLogit ~ pred_challenger$Temp, type = "l", col = "steelblue",
lwd = 2, ylim = c(0,1), xlab = "Temperatura", ylab = "pi.hat", main = "Logit")
plot(pred_challenger$PredPregibin ~ pred_challenger$Temp, type = "l", col = "salmon",
lwd = 2, ylim = c(0,1), xlab = "Temperatura", ylab = "pi.hat", main = "Pregibin")
# Primeiras linhas da base
head(pred_challenger)## Temp PredLogit PredProbit PredClogLog PredPregibin
## 1 31 0.8177744 0.6963991 0.9726399 0.8177731
## 2 32 0.7999113 0.6767297 0.9600989 0.7999099
## 3 33 0.7807665 0.6565534 0.9440667 0.7807651
## 4 34 0.7603387 0.6359209 0.9243209 0.7603373
## 5 35 0.7386447 0.6148865 0.9007993 0.7386432
## 6 36 0.7157212 0.5935085 0.8736052 0.7157197A função Pregibin ajustou probabilidades muito similares a função Logito, mostrando diferença apenas nas casas decimais.
Questão 4
a)
dados <- read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/picloram.csv")
ajuste_picloram <- glm(cbind(kill, total-kill) ~ picloram, family = binomial(link = "logit"), data = dados)
coef(ajuste_picloram)## (Intercept) picloram
## -3.459958 2.656715Assim,
\[-3,46 + 2,66*picloram = log \left(\frac{\hat\pi}{1-\hat\pi}\right)\]
b)
dados$prop <- dados$kill / dados$total
dados$predLogit <- predict(ajuste_picloram, type = "response")
plot(dados$predLogit ~ dados$prop,
xlab = "pi", ylab = "pi.hat")
abline(a = 0, b = 1, col = "red", lty = 2)
ajuste_picloram_pregibin <- glmx(formula = cbind(kill, total-kill) ~ picloram, data = dados,
family = pregibin, xstart = c(0,0), xlink = "logit")
dados$predPregibin <- predict(ajuste_picloram_pregibin, type = "response")
kable(dados)| picloram | rep | kill | total | prop | predLogit | predPregibin |
|---|---|---|---|---|---|---|
| 0.0 | 1 | 0 | 18 | 0.0000000 | 0.0304733 | 0.0000000 |
| 1.1 | 1 | 13 | 24 | 0.5416667 | 0.3687528 | 0.4285714 |
| 2.2 | 1 | 28 | 29 | 0.9655172 | 0.9156620 | 0.8860760 |
| 4.5 | 1 | 24 | 24 | 1.0000000 | 0.9997956 | 1.0000000 |
| 0.0 | 2 | 0 | 29 | 0.0000000 | 0.0304733 | 0.0000000 |
| 1.1 | 2 | 11 | 21 | 0.5238095 | 0.3687528 | 0.4285714 |
| 2.2 | 2 | 18 | 24 | 0.7500000 | 0.9156620 | 0.8860760 |
| 4.5 | 2 | 31 | 31 | 1.0000000 | 0.9997956 | 1.0000000 |
| 0.0 | 3 | 0 | 28 | 0.0000000 | 0.0304733 | 0.0000000 |
| 1.1 | 3 | 9 | 32 | 0.2812500 | 0.3687528 | 0.4285714 |
| 2.2 | 3 | 24 | 26 | 0.9230769 | 0.9156620 | 0.8860760 |
| 4.5 | 3 | 27 | 27 | 1.0000000 | 0.9997956 | 1.0000000 |
É possível identificar um padrão linear satisfatório entre as proporções de ervas daninhas mortas reais e estimadas (função de ligação logito), apesar do nível 1.1 do picloram apresentar um desvio mais elevado.
c)
picloram_0.9 <- (log(0.9 / (1 - 0.9)) - (-3.459958)) / 2.656715
# Estimativa:
picloram_0.9## [1] 2.12939dados$picloram_estimado <- (log(dados$prop / (1 - dados$prop)) - (-3.459958)) / 2.656715
plot(dados$predLogit ~ dados$prop,
xlab = "pi", ylab = "pi.hat")
abline(a = 0, b = 1, col = "red", lty = 2)
# Picloram estimado = pontos azuis
par(new = TRUE)
plot(dados$picloram_estimado ~ dados$prop, col = "blue", yaxt = "n", xaxt = "n", xlab = "", ylab = "", pch = 15)
axis(4)
O picloram estimado parece crescer conforme a proporção de ervas daninhas mortas cresce, de maneira mais acentuada para proporções mais altas.