Questão 1

Mostre a distribuição Poisson Truncada, com função de probabilidade f(y;λ)=λ λy (eλ−1)y

!, em que y=0,1,2,…λ>0 pertence á família exponencial f(y)=exp[ϕ{yθ−b(θ)}+c(y,ϕ)].

exp{log( λy/ (ey−1)y!)}

exp=log()

= 1[ylog(λ)−log(ey−1)]+logy!

-> ϕ=1, θ=log(), b(θ)=log(ey−1), c(y,)=log(y!)$$

Questão 2

Na tabela abaixo descrevemos o número de bactérias sobreviventes em amostras de um produto alimentício segundo o tempo (em minutos) de exposição do produto a uma temperatura de 300◦F

|Número| 175 108 95 82 71 50 49 31 28 17 16 11|

|Tempo| 1 2 3 4 5 6 7 8 9 10 11 12|

Carregando os dados

numero=c(175,   108,    95, 82, 71, 50, 49, 31, 28, 17, 16, 11)
tempo=seq(1:12)
q2 <- data.frame(tempo,numero)
colnames(q2) <- c("Tempo", "Número")
attach(q2)

Histograma Usualmente analisamos o histograma da variável dependente para verificar a melhor distribuição para o ajuste do modelo. Abaixo podemos verificar que a distribuição Poisson aparenta ser adequado para os dados em questão.

hist(Número)

Ajuste

fit1<-glm(Número~Tempo, family = poisson(link=identity))
fit2<-glm(Número~Tempo, family = poisson(link=log))
fit3<-glm(Número~Tempo, family = poisson(link=sqrt))

Análise de Desvio Vamos verificar a adequação dos modelos. Se , o modelo em investigação é aceito.

phi1<-summary(fit1)$dispersion
desvio1<-summary(fit1)$deviance/phi1
q.quadr1<-qchisq(0.95,summary(fit1)$df.residual)

phi2<-summary(fit2)$dispersion
desvio2<-summary(fit2)$deviance/phi2
q.quadr2<-qchisq(0.95,summary(fit2)$df.residual)

phi3<-summary(fit3)$dispersion
desvio3<-summary(fit3)$deviance/phi3
q.quadr3<-qchisq(0.95,summary(fit3)$df.residual)

Pelo comparativo dos desvios o modelo 1 é o único que apresenta resultados não interresante para o ajuste dos dados. Continuaremos nossas analises com o modelo 2 e 3.

library(ggplot2) 
## Warning: package 'ggplot2' was built under R version 4.2.3
library("gridExtra")
va2<-fitted(fit2)
va3<-fitted(fit3) 
y = Número

g2 <- ggplot(q2 , aes(x=va2, y=y)) + ylim(0,200) + xlim(0,200) +
geom_point(size=4, color='#e7ad52') + geom_abline(intercept = 0, slope =
1, color='#001132') + theme(axis.title.x = element_text(colour =
'#001132'), axis.text.x = element_text(colour = '#001132'), axis.title.y
= element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132'))

g3 <- ggplot(q2, aes(x=va3, y=y)) + ylim(0,200) + xlim(0,200) +
geom_point(size=4, color='#e7ad52') + geom_abline(intercept = 0, slope =
1, color='#001132') + theme(axis.title.x = element_text(colour =
'#001132'), axis.text.x = element_text(colour = '#001132'), axis.title.y
= element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132'))

grid.arrange(g2, g3, ncol=2, nrow=1)

Pela análise gráfica o modelo 2, aparenta estar melhor adequada que o modelo 3. Entretanto, o modelo 3 está muito próximo será necessário uma análise de variância para retirar possíveis duvidas.

library(boot)
library("gridExtra")
library(ggplot2)
rd2<-glm.diag(fit2)$rd 
rd3<-glm.diag(fit3)$rd

qq2 <- ggplot(q2, aes(sample=rd2)) + stat_qq(size=4,
color="#e7ad52") + ylim(-4, 4) + xlim(-4, 4) + stat_qq_line(colour =
'#001132') + theme(axis.title.x = element_text(colour = '#001132'),
axis.text.x = element_text(colour = '#001132'), axis.title.y =
element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132'))

qq3 <- ggplot(q2, aes(sample=rd3)) + stat_qq(size=4,
color="#e7ad52") + ylim(-4, 4) + xlim(-4, 4) + stat_qq_line(colour =
'#001132')+ theme(axis.title.x = element_text(colour = '#001132'),
axis.text.x = element_text(colour = '#001132'), axis.title.y =
element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132'))

grid.arrange(qq2, qq3, ncol=2, nrow=1)

Ao analisarmos a função de variância, temos a comprovação de que o modelo 2 estar com a reta muito melhor ajustada que o modelo 3. Por conta disso, usaremos o segundo modelo.

library(knitr) 
ajuste<-summary(fit2)
kable(ajuste$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.3055715 0.0634827 83.57504 0
Tempo -0.2288956 0.0126999 -18.02336 0

\[log(\mu_i) = 5.3 - 0.228*tempo_i\] Ou seja a cada minuto aumentado teremos uma diminuição de 0.228 no log da quantidade média de bactérias.

Análise de Diagnósticos

-Teste de Normalidade

library(nortest)
rd1<-glm.diag(fit2)$rd 
teste1<-lillie.test(rd1)
teste2<-shapiro.test(rd1)

Ao observar os resíduos, temos a informação que são normais, pois seu p-valor dos testes são maiores que 5%, aceitamos a hipótese de que os resíduos apresentam normalidade. Também podemos verificar a normalidade através do gráfico qq-plot.

library(ggplot2) 
qq1 <- ggplot(q2, aes(sample=rd1))+
stat_qq(size=4, color="#e7ad52") + ylim(-4, 4) + xlim(-4, 4) +
stat_qq_line(colour = '#001132') + theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132'))
qq1

  • Diagnósticos Aberrantes
rd1<-glm.diag(fit2)$rd
indice=as.vector(1:length(Número)) 

rd <- ggplot(q2, aes(x=indice,
y=rd1)) + geom_point(size=4, color='#e7ad52') + geom_abline(intercept =
-2, slope = 0, color='#001132')+ geom_abline(intercept = 2, slope = 0,
color='#001132') + ylim(-4,4) + theme(axis.title.x = element_text(colour
= '#001132'), axis.text.x = element_text(colour = '#001132'),
axis.title.y = element_text(colour = '#001132'), axis.text.y =
element_text(colour = '#001132')) 

rd

Como observado, temos apenas um único ponto próximo aos limites.

  • Alavanca
h<-glm.diag(fit2)$h 
indice =as.vector(1:length(Número)) 
h3 <-ggplot(q2, aes(x=indice, y=h)) + geom_point(size=4,
color='#e7ad52') + geom_abline(intercept = 2*(2/length(Número)), slope
= 0, color='#001132') + theme(axis.title.x = element_text(colour =
'#001132'), axis.text.x = element_text(colour = '#001132'), axis.title.y
= element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132')) 

h3

O modelo tem um único ponto acima do limite, podendo ser um possível ponto de alavanca.

  • Influentes
DC<-glm.diag(fit2)$cook

indice=as.vector(1:length(Número))

dc <- ggplot(q2, aes(x=indice, y=DC)) + geom_point(size=4,
color='#e7ad52') + geom_abline(intercept = qchisq(0.1,2)/2, slope = 0,
color='#001132') + theme(axis.title.x = element_text(colour =
'#001132'), axis.text.x = element_text(colour = '#001132'), axis.title.y
= element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132'))

dc

O modelo tem 3 pontos acima dos limites, onde 2 estão muito Destoante dos demais.

#Questão 3

Carregando o banco de dados e nomeando as colunas para facilitar a analise.

library(readr)
questao3  <-read_csv("E:/OneDrive/Documentos/defects.txt", col_names = FALSE)
## Rows: 30 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (4): X1, X2, X3, X4
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
colnames(questao3) <- c("Temperatura", "Densidade",
"Taxa", "Defeitos") 
attach(questao3)

Iremos utilizar a distribuição normal:

fit1<-glm(Defeitos~Temperatura+Densidade+Taxa, family =
gaussian(link=identity))

fit2<-glm(Defeitos~Temperatura+Densidade+Taxa, family =
gaussian(link=log)) 

fit3<-glm(Defeitos~Temperatura+Densidade+Taxa,
family = gaussian(link=inverse))
  • Análise de Desvio

Vamos verificar a adequação dos modelos.

phi1<-summary(fit1)$dispersion
desvio1<-summary(fit1)$deviance/phi1
q.quadr1<-qchisq(0.95,summary(fit1)$df.residual)


phi2<-summary(fit2)$dispersion 
desvio2 <-summary(fit2)$deviance/phi2
q.quadr2<-qchisq(0.95,summary(fit2)$df.residual)


phi3<-summary(fit3)$dispersion
desvio3<-summary(fit3)$deviance/phi3
q.quadr3<-qchisq(0.95,summary(fit3)$df.residual)

Podemos observar que todos os modelos são adequados

library(ggplot2) 
library("gridExtra")

va1<-fitted(fit1) 
va2<-fitted(fit2)
va3<-fitted(fit3) 
y<-Defeitos


g1 <- ggplot(questao3 , aes(x=va1, y=y)) +
geom_point(size=4, color='#e7ad52') + geom_abline(intercept = 0, slope =
1, color='#001132') + ggtitle("Modelo 1") + theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132')) 

g2<- ggplot(questao3, aes(x=va2,
y=y)) + geom_point(size=4, color='#e7ad52') + geom_abline(intercept = 0,
slope = 1, color='#001132') + ggtitle("Modelo 2") + theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132')) 

g3 <- ggplot(questao3, aes(x=va3,
y=y)) + geom_point(size=4, color='#e7ad52') + geom_abline(intercept = 0,
slope = 1, color='#001132') + ggtitle("Modelo 3") + theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132'))

grid.arrange(g1, g2, g3, ncol=3, nrow=1)

Observamos que o modelo 1 e o modelo 2 são os melhores ajustes diante da reta. Entretanto, entre os dois o modelo um aparenta ter melhores resultados.

Escolha da função de variância Gráfico de versus valores ajustados.

library("gridExtra")
library(ggplot2) 
library(boot) 

rd1<-glm.diag(fit1)$rd 
rd2<-glm.diag(fit2)$rd
rd3<-glm.diag(fit3)$rd

q1 <- ggplot(questao3, aes(sample=rd1))+
ggtitle("Modelo 1") + stat_qq(size=4, color="#e7ad52") + ylim(-4, 4) +
xlim(-4, 4) + stat_qq_line(colour = '#001132')+ theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132')) 

q2 <- ggplot(questao3,
aes(sample=rd2)) + ggtitle("Modelo 2") + stat_qq(size=4,
color="#e7ad52") + ylim(-4, 4) + xlim(-4, 4) + stat_qq_line(colour =
'#001132')+ theme(axis.title.x = element_text(colour = '#001132'),
axis.text.x = element_text(colour = '#001132'), axis.title.y =
element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132')) 

q3 <- ggplot(questao3, aes(sample=rd3)) + stat_qq(size=4,
color="#e7ad52") + ylim(-4, 4) + xlim(-4, 4) + ggtitle("Modelo 3") +
stat_qq_line(colour = '#001132')+ theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132')) 

grid.arrange(q1, q2, q3, ncol=3, nrow=1)

Por fim, através do gráfico acima, fica mais evidente a escolha do modelo 1 como o mais apropriado para nossos dados.

Seleção de variáveis Com a função identidade sendo escolha mais apropriada, iremos rodar o algoritmo de stepwise para tentarmos identificar as variáveis mais significantes para nosso modelo.

menor<-glm(Defeitos~1, family = gaussian(link=identity))

maior<-glm(Defeitos~Temperatura+Densidade+Taxa,family =
gaussian(link=identity)) 
step(menor, scope = list(upper=maior),
direction="both") 
## Start:  AIC=266.06
## Defeitos ~ 1
## 
##               Df Deviance    AIC
## + Temperatura  1   1496.9 208.43
## + Densidade    1   1612.5 210.67
## + Taxa         1   2364.3 222.15
## <none>            10923.9 266.06
## 
## Step:  AIC=208.43
## Defeitos ~ Temperatura
## 
##               Df Deviance    AIC
## + Densidade    1   1354.8 207.44
## + Taxa         1   1389.7 208.21
## <none>             1496.9 208.43
## - Temperatura  1  10923.9 266.06
## 
## Step:  AIC=207.44
## Defeitos ~ Temperatura + Densidade
## 
##               Df Deviance    AIC
## <none>             1354.8 207.44
## - Densidade    1   1496.8 208.43
## + Taxa         1   1314.4 208.53
## - Temperatura  1   1612.5 210.67
## 
## Call:  glm(formula = Defeitos ~ Temperatura + Densidade, family = gaussian(link = identity))
## 
## Coefficients:
## (Intercept)  Temperatura    Densidade  
##      46.238       18.050       -2.327  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       10920 
## Residual Deviance: 1355  AIC: 207.4

Modelo Final Através das análises de desvio, análises gráficas e seleção de variáveis foi escolhido o modelo 1 que faz uso da funçao identidade e as seguintes variáveis:

library(knitr) 
fit1<-glm(Defeitos ~ Temperatura + Densidade, family =
gaussian(link=identity)) 
ajuste<-summary(fit1)
kable(ajuste$coefficients)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 46.237985 52.065528 0.8880729 0.3823449
Temperatura 18.049940 7.965269 2.2660803 0.0316738
Densidade -2.327472 1.383311 -1.6825371 0.1039958

Isto é, fixando a densidade, quando quando acrescido uma unidade de temperatura teremos um aumento no número médio de defeitos. E quando fixado a temperatura, quando acrescido uma unidade de densidade teremos uma diminuição no número médio de defeitos.

-Análise de Diagnósticos -Teste de Normalidade

library(nortest)
rd1<-glm.diag(fit1)$rd 
teste1<-lillie.test(rd1)
teste2<-shapiro.test(rd1)

Os resíduos são normais,como o p-valor dos testes aplicados são maiores que 5%.

library(ggplot2) 

q1 <- ggplot(questao3, aes(sample=rd1))+
stat_qq(size=4, color="#e7ad52") + ylim(-4, 4) + xlim(-4, 4) +
stat_qq_line(colour = '#001132') + theme(axis.title.x =
element_text(colour = '#001132'), axis.text.x = element_text(colour =
'#001132'), axis.title.y = element_text(colour = '#001132'), axis.text.y
= element_text(colour = '#001132'))
q1

Diagnósticos Aberrantes

rd1<-glm.diag(fit1)$rd

indice=as.vector(1:length(Defeitos)) 
rd <- ggplot(questao3,
aes(x=indice, y=rd1)) + geom_point(size=4, color='#e7ad52') +
geom_abline(intercept = -2, slope = 0, color='#001132')+
geom_abline(intercept = 2, slope = 0, color='#001132') + ylim(-4,4) +
theme(axis.title.x = element_text(colour = '#001132'), axis.text.x =
element_text(colour = '#001132'), axis.title.y = element_text(colour =
'#001132'), axis.text.y = element_text(colour = '#001132')) 
rd

O modelo apresentou apenas um único ponto acima dos limites e um único ponto abaixo dos limites.

  • Alavanca
h<-glm.diag(fit1)$h

indice=as.vector(1:length(Defeitos)) 
h3<- ggplot(questao3, aes(x=indice, y=h)) + geom_point(size=4,
color='#e7ad52') + geom_abline(intercept =
2*(length(coef(fit1))/length(Defeitos)), slope = 0, color='#001132') +
theme(axis.title.x = element_text(colour = '#001132'), axis.text.x =
element_text(colour = '#001132'), axis.title.y = element_text(colour =
'#001132'), axis.text.y = element_text(colour = '#001132')) 
h3

Para o gráfico de pontos de alavanca, o modelo apresentou apenas um único ponto acima do limite estabelecido.

-Influentes

DC<-glm.diag(fit1)$cook

indice=as.vector(1:length(Defeitos)) 
dc <- ggplot(questao3,
aes(x=indice, y=DC)) + geom_point(size=4, color='#e7ad52') +
geom_abline(intercept = qchisq(0.05,length(coef(fit1)))/2, slope = 0,
color='#001132') + theme(axis.title.x = element_text(colour =
'#001132'), axis.text.x = element_text(colour = '#001132'), axis.title.y
= element_text(colour = '#001132'), axis.text.y = element_text(colour =
'#001132'))

dc

O modelo tem 2 pontos acima dos limites, talvez sendo necessária uma investigação mais aprofundada.