a)

Stormchaser = read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/Stormchaser.csv")

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 = 'Tornadoes', ylab = 'Prob', 
     lwd = 2, xaxt = 'n', ylim = c(0,1), col = 'steelblue4')
axis(side = 1, at = 0:8)
lines(x = y + 0.1, y = rel.freq2, type = 'h',
      lwd = 2, lty = 'solid', col = 'indianred')
grid()

legend("topright", legend = c("Poisson", "Valores Observados"), 
       col = c("steelblue4", "indianred"), pch = c(1, 16))

Com base no gráfico apresentado, é possível notar que os valores observados da variável “Tornadoes” exibem uma frequência significativa de ocorrências igual a zero. Esse padrão sugere a possibilidade de explorar um modelo com suporte para inflação de zeros, o que pode ser considerado como uma alternativa.

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)

O parâmetro buscaria quantificar a probabilidade de não observarmos nenhum tornados, mesmo que as condições climaticas estiverem adequadas para um.

d)

Abaixo ajustaremos os dois modelos.

require(pscl)
#Ajustando os modelos
ajuste_zip <- zeroinfl(tornadoes~1|1, dist = 'poisson', data = Stormchaser)
ajuste_hurdle <- hurdle(tornadoes~1, dist = 'poisson', data=Stormchaser)

Agora iremos calcular as estimativas e IC de \(\hat\mu\) e \(\hat\pi\) do ajuste ZIP.

#ZIP mu_hat
round(coef(ajuste_zip)[1],2); round(confint(ajuste_zip)[1,],2)

## count_(Intercept) 
##              0.41

##  2.5 % 97.5 % 
##  -0.07   0.89

#ZIP Pi_hat
round(coef(ajuste_zip)[2],2); round(confint(ajuste_zip)[2,],2)

## zero_(Intercept) 
##             0.95

##  2.5 % 97.5 % 
##   0.26   1.63

Agora faremos a mesma coisa para o ajuste Hurdle.

#Hurdle mu_hat
round(coef(ajuste_hurdle)[1],2); round(confint(ajuste_hurdle)[1,],2)

## count_(Intercept) 
##              0.41

##  2.5 % 97.5 % 
##  -0.07   0.89

#Hurdle Pi_hat
round(coef(ajuste_hurdle)[2],2); round(confint(ajuste_hurdle)[2,],2)

## zero_(Intercept) 
##            -1.28

##  2.5 % 97.5 % 
##  -1.85  -0.71

Podemos observar que para ambos os modelos ajustados, o IC para \(\hat\mu\) inclui o valor zero, sendo um indicativo de que o parâmetro pode ser não significativo.

a)

O número de clientes que esperam na fila na loja do Starbucks. Com base na amostra coletada na loja de Lincoln, NE, entre às 08:00 e 08:30 da manhã

b)

starbucks <-  read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/starbucks.csv" )
starbucks$Day <- factor(starbucks$Day, 
                       levels = c('Monday','Tuesday','Wednesday','Thursday','Friday'))
boxplot(starbucks$Count~starbucks$Day, col = 'indianred4', xlab= "Dia da Semana", ylab = "Número de Clientes",
        main = "Boxplot Clientes por Dia da Semana")
points(starbucks$Count, col = 'royalblue1')

Com base no gráfico apresentado, não parece plausível assumir que o número médio de clientes é igual todos os dias da semana. Segunda-feira (Monday) aparece como o dia mais discrepante dentre os demais.

c)

Abaixo ajustamos um modelo com um termo offset que considera o número de vezes em que cada dia da semana ocorreu na base.

#Definindo o termo offset.
starbucks$t <- as.numeric(table(starbucks$Day))
ajuste <- glm(Count~Day, family = poisson(link = 'log'), offset = log(t),data = starbucks)
summary(ajuste)

## 
## Call:
## glm(formula = Count ~ Day, family = poisson(link = "log"), data = starbucks, 
##     offset = log(t))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4641  -0.8832  -0.5099   0.9392   3.2908  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.1394     0.3536  -3.223 0.001269 ** 
## DayTuesday     0.9651     0.4155   2.323 0.020188 *  
## DayWednesday   1.0986     0.4082   2.691 0.007123 ** 
## DayThursday    1.3218     0.3979   3.322 0.000895 ***
## DayFriday      1.1394     0.4062   2.805 0.005030 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 87.432  on 24  degrees of freedom
## Residual deviance: 72.431  on 20  degrees of freedom
## AIC: 150.9
## 
## Number of Fisher Scoring iterations: 5

Realizando o LRT

require(car)
Anova(ajuste, 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

Com base nesse resultado, podemos observar que o dia há evidências significativas, a um um nível de confiança de 95%, que o dia da semana é significativo.