# Questão 1
Stormchaser = read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/Stormchaser.csv")
# a)
rel.freq <- table(Stormchaser$tornadoes) / length(Stormchaser$tornadoes)
rel.freq2 <- c(rel.freq, rep(0, times = 2))
x.bar <- mean(Stormchaser$tornadoes) ; y <- 0:7
prob <- round(dpois(x=y, lambda = x.bar), 4)
data.frame(y, prob, rel.freq = rel.freq2)
## y prob rel.freq
## 1 0 0.6569 0.78260870
## 2 1 0.2761 0.13043478
## 3 2 0.0580 0.04347826
## 4 3 0.0081 0.01449275
## 5 4 0.0009 0.01449275
## 6 5 0.0001 0.01449275
## 7 6 0.0000 0.00000000
## 8 7 0.0000 0.00000000
plot(x = y - 0.1, y = prob, type = 'h',
xlab = 'Tornados', ylab = 'Probabilidade',
lwd = 2, xaxt = 'n', ylim = c(0,1), col = 'steelblue')
axis(side = 1, at = 0:8)
lines(x = y + 0.1, y = rel.freq2, type = 'h',
lwd = 2, lty = 'solid', col = '#ff5454')
grid()
# Linha azul = observado
# Linha rosa = Poisson
# b)
# O modelo Hurdle pressupõe dois processos distintos,
# onde um determina a presença ou ausência de zeros e o
# outro gera valores positivos. O ZIP assume que zeros
# podem surgir de duas maneiras dentro do mesmo ajuste
# (zeros estruturais e zeros aleatórios).
# c)
# Nesse contexto, o parâmetro quantificaria a probabilidade
# de observarmos zero tornados mesmo quando as condições climáticas
# favorecem a ocorrência dos mesmos.
# d)
library(pscl)
zip <- zeroinfl(tornadoes ~ 1 | 1, dist = "poisson", data = Stormchaser)
hurdle <- hurdle(tornadoes ~ 1, dist = "poisson", data = Stormchaser)
# Estimativa e intervalo de confiança mu.hat ZIP
round(coef(zip)[1], 2) ; round(confint(zip)[1,], 2)
## count_(Intercept)
## 0.41
## 2.5 % 97.5 %
## -0.07 0.89
# Estimativa e intervalo de confiança mu.hat Hurdle
round(coef(hurdle)[1], 2) ; round(confint(hurdle)[1,], 2)
## count_(Intercept)
## 0.41
## 2.5 % 97.5 %
## -0.07 0.89
# Estimativa e intervalo de confiança exp(pi.hat) ZIP
round(exp(coef(zip)[2]), 2) ; round(exp(confint(zip)[2,]), 2)
## zero_(Intercept)
## 2.58
## 2.5 % 97.5 %
## 1.30 5.11
# Estimativa e intervalo de confiança exp(pi.hat) Hurdle
round(exp(coef(hurdle)[2]), 2) ;round(exp(confint(hurdle)[2,]), 2)
## zero_(Intercept)
## 0.28
## 2.5 % 97.5 %
## 0.16 0.49
# Questão 2
starbucks = read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/starbucks.csv")
# a)
# População de inferência: Filiais da starbucks.
# Amostra: Filial de Lincoln, NE, entre às 08:00 e 08:30 da manhã.
# b)
starbucks$Day <- factor(starbucks$Day,
levels = c('Monday','Tuesday','Wednesday','Thursday','Friday'))
plot(Count ~ as.numeric(starbucks$Day), data = starbucks,
xlab = 'Dia da semana', ylab = 'Clientes na fila',
col = '#ff5454', lwd = 1.5, xaxt = 'n')
axis(1, at = 1:5, labels = c('Segunda', 'Terça', 'Quarta', 'Quinta', 'Sexta'))
# Não parece plausível que a verdadeira contagem média seja constante ao longo dos dias,
# pois é visível um padrão crescente ao decorrer da semana, apesar dos pontos dos mesmos dias
# nem sempre apresentarem o mesmo comportamento (principalmente na sexta).
# c)
starbucks$t <- 5 # Total de dias
ajuste1 <- glm(Count ~ Day, family = poisson(link = "log"), offset = log(t), data = starbucks)
# Estimação do modelo
round(coef(ajuste1), 2)
## (Intercept) DayTuesday DayWednesday DayThursday DayFriday
## -1.14 0.97 1.10 1.32 1.14
library(car)
Anova(ajuste1, test = 'LR')
## Analysis of Deviance Table (Type II tests)
##
## Response: Count
## LR Chisq Df Pr(>Chisq)
## Day 15.002 4 0.004698 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Dia da semana é significativo
# Questão 3
BladderCancer = read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/BladderCancer.csv")
ajuste_b_c <- glm(Tumors ~ Size, data = BladderCancer, family = poisson(link = "log"), offset = log(Time))
summary(ajuste_b_c)
##
## Call:
## glm(formula = Tumors ~ Size, family = poisson(link = "log"),
## data = BladderCancer, offset = log(Time))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3394 0.1768 -13.234 <2e-16 ***
## Size 0.3080 0.3062 1.006 0.315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 29.850 on 30 degrees of freedom
## Residual deviance: 28.876 on 29 degrees of freedom
## AIC: 103.69
##
## Number of Fisher Scoring iterations: 5
# O tamanho do tumor primário inicial não está relacionado ao número de tumores (não é significativo no modelo)
# Questão 4
library(glmnet)
dt <- read.csv("http://leg.ufpr.br/~lucambio/ADC/dt.csv")
xx <- as.matrix(dt[, c("hosp", "numchron", "gender", "school", "health_excellent", "health_poor")])
yy <- dt$privins
lasso.fit <- glmnet (y = yy , x = xx , family = "binomial")
las.lambda <- lasso.fit$lambda
length(las.lambda)
## [1] 52
las.coefs <- Matrix(rbind(as.matrix(coef(lasso.fit)), las.lambda), sparse = TRUE)
round(las.coefs[,c(1:5)], digits = 3)
## 8 x 5 sparse Matrix of class "dgCMatrix"
## s0 s1 s2 s3 s4
## (Intercept) 1.245 1.053 0.883 0.731 0.593
## hosp . . . . .
## numchron . . . . .
## gender . . . . .
## school . 0.019 0.036 0.051 0.065
## health_excellent . . . . .
## health_poor . . . . .
## las.lambda 0.138 0.126 0.115 0.104 0.095
set.seed (123)
cv.lasso.fit <- cv.glmnet (y = yy , x = xx , family = "binomial")
plot(cv.lasso.fit)
grid()
predict.las <- predict (cv.lasso.fit , newx = xx , type = "response")
head(cbind(dt$privins,round(predict.las,digits = 3)))
## lambda.1se
## [1,] 1 0.659
## [2,] 1 0.788
## [3,] 0 0.766
## [4,] 1 0.511
## [5,] 1 0.659
## [6,] 0 0.668