rm(list=ls(all=T))
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 561883 30.1 1283252 68.6 644254 34.5
Vcells 997387 7.7 8388608 64.0 1635074 12.5Primeira Prova - Regressão 2
Aluno: Paulo Manoel da Silva Junior
Matrícula: 20190041314
Questão 1
- Mostre a distribuição Poisson Truncada, com função de probabilidade
\[f(y;\lambda) = \frac{\lambda^{y}}{(e^\lambda - 1)y!},\]
em que \(y = 0,1,2,3,...., \quad \lambda >0\) pertence à família exponencial \((f(y) = \exp[\phi\left \{ y\theta - b(\theta) \right \} + c(y;\theta)])\).
Aplicando a família exponencial e as propriedades do log, temos:
\[= \exp \left \{ y \log \lambda - \log(\exp^\lambda-1) + \log y!\right \}\]
\[= \ \exp \left\{ 1 \times ( y \log \lambda - \log (\exp^\lambda - 1) + \log y! \right \}\]
Logo, vemos que: \(\phi = 1\), \(\theta = \log \lambda\), \(b(\theta) = \log (\exp^\lambda - 1)\) e \(c(y, \phi) = \log y!\).
Sendo assim, conclui-se que a Poisson Truncada pertence à família exponencial e pode ser utilizado no MLG.
Questão 2
- Descrevemos na tabela abaixo o
númerode bactérias sobreviventes em amostras de um produto alimentício segundo otempo(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 |
numero <- c(175,108,95,82,71,50,49,31,28,17,16,11)
tempo <- c(1:12)
data <- data.frame(numero, tempo)Estatística Descritiva
skimr::skim(data)| Name | data |
| Number of rows | 12 |
| Number of columns | 2 |
| _______________________ | |
| Column type frequency: | |
| numeric | 2 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| numero | 0 | 1 | 61.08 | 48.26 | 11 | 25.25 | 49.5 | 85.25 | 175 | ▇▅▅▁▂ |
| tempo | 0 | 1 | 6.50 | 3.61 | 1 | 3.75 | 6.5 | 9.25 | 12 | ▇▅▅▅▇ |
DescTools::Desc(data, plotit = F, digits = 3)------------------------------------------------------------------------------
Describe data (data.frame):
data frame: 12 obs. of 2 variables
12 complete cases (100.0%)
Nr ColName Class NAs Levels
1 numero numeric .
2 tempo integer .
------------------------------------------------------------------------------
1 - numero (numeric)
length n NAs unique 0s mean meanCI'
12 12 0 = n 0 61.083 30.422
100.0% 0.0% 0.0% 91.745
.05 .10 .25 median .75 .90 .95
13.750 16.100 25.250 49.500 85.250 106.700 138.150
range sd vcoef mad IQR skew kurt
164.000 48.258 0.790 48.184 60.000 0.935 -0.065
value freq perc cumfreq cumperc
1 11 1 8.3% 1 8.3%
2 16 1 8.3% 2 16.7%
3 17 1 8.3% 3 25.0%
4 28 1 8.3% 4 33.3%
5 31 1 8.3% 5 41.7%
6 49 1 8.3% 6 50.0%
7 50 1 8.3% 7 58.3%
8 71 1 8.3% 8 66.7%
9 82 1 8.3% 9 75.0%
10 95 1 8.3% 10 83.3%
11 108 1 8.3% 11 91.7%
12 175 1 8.3% 12 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
2 - tempo (integer)
length n NAs unique 0s mean meanCI'
12 12 0 = n 0 6.500 4.209
100.0% 0.0% 0.0% 8.791
.05 .10 .25 median .75 .90 .95
1.550 2.100 3.750 6.500 9.250 10.900 11.450
range sd vcoef mad IQR skew kurt
11.000 3.606 0.555 4.448 5.500 0.000 -1.502
value freq perc cumfreq cumperc
1 1 1 8.3% 1 8.3%
2 2 1 8.3% 2 16.7%
3 3 1 8.3% 3 25.0%
4 4 1 8.3% 4 33.3%
5 5 1 8.3% 5 41.7%
6 6 1 8.3% 6 50.0%
7 7 1 8.3% 7 58.3%
8 8 1 8.3% 8 66.7%
9 9 1 8.3% 9 75.0%
10 10 1 8.3% 10 83.3%
11 11 1 8.3% 11 91.7%
12 12 1 8.3% 12 100.0%
' 95%-CI (classic)Histograma
É observardo que a variável resposta é o número de bactérias sobreviventes.
ggplot2::ggplot(data, ggplot2::aes(x = numero))+
ggplot2::geom_histogram(fill='white',
color = "blue",
breaks = hist(data$numero, plot = F)$breaks) + ggplot2::xlab("Número") + ggplot2::ylab("Frequência") + ggplot2::ggtitle("Histograma do Número de Bactérias")
hist(data$numero, xlab = "Número", ylab = "Frequência", main = "Histograma do número de bactérias")
Ajuste dos modelos com as funções de ligação
- Realize o ajuste da poisson com as possíveis funções de ligação e decida, entre elas, qual a função de ligação é melhor. Justifique sua escolha.
fit1<-glm(numero ~ tempo, family = poisson(link=log), data = data)
fit2<-glm(numero ~ tempo, family = poisson(link=sqrt), data = data)
fit3<-glm(numero ~ tempo, family = poisson(link=identity), data = data)Análise do Desvio
Vamos verificar a adequação dos modelos. Se \(\frac{D (y;\mu)}{\phi} \leq \chi^2_{n-p;1-\alpha}\) , o modelo em investigação é aceito.
Fixando o \(\alpha = 0.05\)
phi1<-summary(fit1)$dispersion
desvio1<-summary(fit1)$deviance/phi1 # log
q.quadr1<-qchisq(0.95,summary(fit1)$df.residual)
phi2<-summary(fit2)$dispersion
desvio2<-summary(fit2)$deviance/phi2 # sqrt
q.quadr2<-qchisq(0.95,summary(fit2)$df.residual)
phi3<-summary(fit3)$dispersion
desvio3<-summary(fit3)$deviance/phi3 # identidade
q.quadr3<-qchisq(0.95,summary(fit3)$df.residual)Verificando se todos os modelos foram aceitos:
| Modelo 1 | Modelo 2 | Modelo 3 | |
|---|---|---|---|
| Desvio/\(\phi\) | 8.4215055 | 15.6463184 | 45.9830251 |
| \(\chi^2\) | 18.3070381 | 18.3070381 | 18.3070381 |
É observardo que o modelo com a função de ligação identidade não é aceito, pois o resultado do \(D (y;\mu)/ \phi\) foi maior que o quantil, ou seja: 45.983 > 18.307. Logo, prosseguiremos analisando apenas a log e sqrt.
Escolha da função de ligação
Vamos então fazer a escolha da função de ligação utilizando o gráfico dos valores observados \((y_i)\) versus valores estimados \((\hat\mu_i)\).
va1<-fitted(fit1)
va2<-fitted(fit2)
y<-numerog1 <- ggplot2::ggplot(data, ggplot2::aes(x=va1,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,180) + ggplot2::xlim(0,180) + ggplot2::ggtitle("Log")
G1 <- g1+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
g2 <- ggplot2::ggplot(data, ggplot2::aes(x=va2,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,180) + ggplot2::xlim(0,180) +
ggplot2::ggtitle("Sqrt")
G2 <- g2+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))gridExtra::grid.arrange(G1,G2, nrow = 1, ncol = 2)
De acordo com a análise gráfica é observado que a função de ligação Log visualmente parece estar melhor adequada do que com a função sqrt.
Escolha da função de Variância
- Gráfico de \(r^*_D\) versus valores ajustados.
rd1 <- boot::glm.diag(fit1)$rd
rd2 <- boot::glm.diag(fit2)$rdq1 <- ggplot2::ggplot(data, ggplot2::aes(sample = rd1))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Log")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
q2 <- ggplot2::ggplot(data, ggplot2::aes(sample = rd2))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Sqrt")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))gridExtra::grid.arrange(q1, q2, nrow = 1, ncol = 2)
Conforme observado na análise gráfica a função de variância Log leva vantagem no exemplo considerado a função de variância sqrt, onde boa parte dos pontos estão mais próximos da reta do que com a função sqrt.
fit <- glm(numero~tempo, family = poisson(log), data = data)
summary(fit)
Call:
glm(formula = numero ~ tempo, family = poisson(log), data = data)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7703 -0.5715 -0.1019 0.5496 1.2794
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.30557 0.06348 83.58 <2e-16 ***
tempo -0.22890 0.01270 -18.02 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 393.6292 on 11 degrees of freedom
Residual deviance: 8.4215 on 10 degrees of freedom
AIC: 80.182
Number of Fisher Scoring iterations: 4- É escolhido o modelo ajustado com a função de ligação log, e a equação que representa o modelo ajustado ficou como:
\[log_{\mu_i} = 5.3056 -0.2289\times tempo_i\]
Conforme temos o acréscimo de uma unidade na quantidade de tempo, vamos ter uma diminuição de \(0.22890\) valor do log da quantidade média de bactérias.
ajuste<-summary(fit)
knitr::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 |
- Realize uma análise de diagnóstico para o melhor modelo obtido em (a). Analise os resultados obtidos.
Análise Diagóstico
Pontos aberrantes
indice=as.vector(1:length(numero))
rd <- boot::glm.diag(fit)$rd
rd1 <- ggplot2::ggplot(data, ggplot2::aes(x=indice, y=rd)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = -2, slope = 0, color='#001132')+ ggplot2::geom_abline(intercept = 2, slope = 0, color='#001132') + ggplot2::ylim(-4,4) + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
rd1
Visualmente, é observado que um ponto parece estar fora dos limites estabelecidos, sendo assim um possível ponto aberrante, no qual é a observação de número 2.
data[which(rd>2 | rd < -2),] numero tempo
2 108 2Pontos de Alavanca
h<- boot::glm.diag(fit)$h
h3 <- ggplot2::ggplot(data, ggplot2::aes(x=indice, y=h)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(fit))/length(numero)), slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
h3
Foi encontrado um ponto de alavanca, que é a observação de número 1. Sendo as observações deste ponto de alavanca
data[which(h>2*(length(coef(fit))/length(numero))),] numero tempo
1 175 1Pontos de Influência
DC<-boot::glm.diag(fit)$cook
dc <- ggplot2::ggplot(data, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.05,length(coef(fit)))/2, slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
dc
Como o nível de significância foi fixado bem mais acima sendo 5%, encontramos 4 pontos de influência, que são as observações de número 1, 2, 7 e 10. E os valores associados estão descritos a seguir.
data[which(DC>qchisq(0.05,length(coef(fit)))/2),] numero tempo
1 175 1
2 108 2
7 49 7
10 17 10De acordo com a natureza do problema, pode ser averiguado para se tomar a melhor decisão para com estes pontos.
Questão 3
- No banco de dados
defects.txtindica a temperatura do processo de produção, a densidade do produto, taxa de produção e a média do número de defeitos nos produtos (variável resposta).
rm(list=ls(all=T))
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 2279570 121.8 4367068 233.3 4367068 233.3
Vcells 3916767 29.9 8388608 64.0 6614368 50.5setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/PROVA 1")
dados <- read.table("defects.txt", sep = ",", col.names = c("temperatura", "densidade", "taxa", "defeitos"), header = F)Estatística Descritiva
skimr::skim(dados)| Name | dados |
| Number of rows | 30 |
| Number of columns | 4 |
| _______________________ | |
| Column type frequency: | |
| numeric | 4 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| temperatura | 0 | 1 | 2.20 | 0.58 | 0.97 | 1.76 | 2.26 | 2.72 | 3.02 | ▂▅▃▆▇ |
| densidade | 0 | 1 | 25.29 | 3.36 | 19.45 | 22.82 | 24.62 | 27.48 | 32.19 | ▃▇▅▃▃ |
| taxa | 0 | 1 | 236.52 | 26.05 | 177.70 | 221.28 | 238.50 | 253.88 | 281.90 | ▂▃▇▇▃ |
| defeitos | 0 | 1 | 27.14 | 19.41 | 0.10 | 10.62 | 26.15 | 43.85 | 60.80 | ▇▅▃▆▅ |
DescTools::Desc(dados, plotit = F, digits = 3)------------------------------------------------------------------------------
Describe dados (data.frame):
data frame: 30 obs. of 4 variables
30 complete cases (100.0%)
Nr ColName Class NAs Levels
1 temperatura numeric .
2 densidade numeric .
3 taxa numeric .
4 defeitos numeric .
------------------------------------------------------------------------------
1 - temperatura (numeric)
length n NAs unique 0s mean meanCI'
30 30 0 29 0 2.203 1.985
100.0% 0.0% 0.0% 2.421
.05 .10 .25 median .75 .90 .95
1.252 1.459 1.765 2.260 2.720 2.857 2.936
range sd vcoef mad IQR skew kurt
2.050 0.583 0.265 0.719 0.955 -0.416 -1.012
lowest : 0.970, 1.090, 1.450, 1.460, 1.500
highest: 2.840, 2.850, 2.920, 2.950, 3.020
' 95%-CI (classic)
------------------------------------------------------------------------------
2 - densidade (numeric)
length n NAs unique 0s mean meanCI'
30 30 0 = n 0 25.289 24.034
100.0% 0.0% 0.0% 26.543
.05 .10 .25 median .75 .90 .95
20.870 21.536 22.817 24.615 27.475 30.079 31.459
range sd vcoef mad IQR skew kurt
12.740 3.359 0.133 3.113 4.658 0.448 -0.737
lowest : 19.450, 20.650, 21.140, 21.580, 22.300
highest: 29.420, 30.010, 30.700, 32.080, 32.190
' 95%-CI (classic)
------------------------------------------------------------------------------
3 - taxa (numeric)
length n NAs unique 0s mean meanCI'
30 30 0 29 0 236.517 226.789
100.0% 0.0% 0.0% 246.244
.05 .10 .25 median .75 .90 .95
186.530 206.390 221.275 238.500 253.875 266.390 273.040
range sd vcoef mad IQR skew kurt
104.200 26.051 0.110 24.018 32.600 -0.488 -0.402
lowest : 177.700, 181.400, 192.800, 207.900, 210.800
highest: 260.000, 265.700, 272.600, 273.400, 281.900
' 95%-CI (classic)
------------------------------------------------------------------------------
4 - defeitos (numeric)
length n NAs unique 0s mean meanCI'
30 30 0 = n 0 27.143 19.896
100.0% 0.0% 0.0% 34.391
.05 .10 .25 median .75 .90 .95
0.785 2.130 10.625 26.150 43.850 51.350 57.215
range sd vcoef mad IQR skew kurt
60.700 19.408 0.715 25.427 33.225 0.110 -1.403
lowest : 0.100, 0.200, 1.500, 2.200, 2.800
highest: 49.700, 50.900, 55.400, 58.700, 60.800
' 95%-CI (classic)Histograma
ggplot2::ggplot(dados, ggplot2::aes(x = defeitos))+
ggplot2::geom_histogram(fill='white',
color = "blue",
breaks = hist(dados$defeitos, plot = F)$breaks) + ggplot2::xlab("Número médio de defeitos") + ggplot2::ylab("Frequência") + ggplot2::ggtitle("Histograma do Número médio de defeitos")
hist(dados$defeitos, xlab = "Número médio de defeitos", ylab = "Frequência", main = "Histograma do número médio de defeitos")
Ajuste dos modelos com as funções de ligação
- Realize o ajuste da normal com as possíveis funções de ligação e decida, entre elas, qual a função de ligação é melhor. Justifique sua escolha.
attach(dados)
fit1<-glm(defeitos ~ temperatura + densidade + taxa, family = gaussian(link=identity), data = dados)
fit2<-glm(defeitos ~ temperatura + densidade + taxa, family = gaussian(link=log), data = dados)
fit3<-glm(defeitos ~ temperatura + densidade + taxa, family = gaussian(link=inverse), data = dados)Análise do Desvio
Vamos verificar a adequação dos modelos. Se \(\frac{D(y;\mu)}{\phi} \leq \chi^2_{n-p;1-\alpha}\) , o modelo em investigação é aceito.
Fixando o \(\alpha = 0.05\)
phi1<-summary(fit1)$dispersion
desvio1<-summary(fit1)$deviance/phi1 # identidade
q.quadr1<-qchisq(0.95,summary(fit1)$df.residual)
phi2<-summary(fit2)$dispersion
desvio2<-summary(fit2)$deviance/phi2 # log
q.quadr2<-qchisq(0.95,summary(fit2)$df.residual)
phi3<-summary(fit3)$dispersion
desvio3<-summary(fit3)$deviance/phi3 # inversa
q.quadr3<-qchisq(0.95,summary(fit3)$df.residual)Verificando se todos os modelos foram aceitos:
| Modelo 1 | Modelo 2 | Modelo 3 | |
|---|---|---|---|
| Desvio/\(\phi\) | 26 | 26.0006548 | 26.0013959 |
| \(\chi^2\) | 38.8851387 | 38.8851387 | 38.8851387 |
É verificado que todos os modelos foram aceitos
Escolha da função de ligação
Vamos então fazer a escolha da função de ligação utilizando o gráfico dos valores observados \((y_i)\) versus valores estimados \((\hat\mu_i)\).
va1<-fitted(fit1)
va2<-fitted(fit2)
va3<-fitted(fit3)
y<-defeitosg1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va1,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,100) + ggplot2::xlim(-20,80) + ggplot2::ggtitle("Identidade")
G1 <- g1+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
g2 <- ggplot2::ggplot(dados, ggplot2::aes(x=va2,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,100) + ggplot2::xlim(0,80) +
ggplot2::ggtitle("Log")
G2 <- g2+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
g3 <- ggplot2::ggplot(dados, ggplot2::aes(x=va3,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,100) + ggplot2::xlim(0,80)+
ggplot2::ggtitle("Inversa")
G3 <- g3+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))gridExtra::grid.arrange(G1,G2,G3, nrow = 1, ncol = 3)
De acordo com a análise gráfica é observado que a função de identidade visualmente parece estar melhor adequada as outras duas funções de ligação
Escolha da função de Variância
- Gráfico de \(r^*_D\) versus valores ajustados.
rd1 <- boot::glm.diag(fit1)$rd
rd2 <- boot::glm.diag(fit2)$rd
rd3 <- boot::glm.diag(fit3)$rdq1 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd1))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Identidade")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
q2 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd2))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Log")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
q3 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd3))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Inversa")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))gridExtra::grid.arrange(q1, q2, q3, nrow = 1, ncol = 3)
Vamos utilizar a função de ligação identidade no modelo ajustado, pois, conforme a análise gráfica pareceu estar melhor ajustada principalmente levando em consideração o gráfico da função de ligação.
Seleção de Variáveis
Vamos considerar a utilização da regressão StepWise para encontrar o melhor modelo de acordo com o menor AIC.
menor <- glm(defeitos ~ 1, family = gaussian(link = identity), data = dados)
maior <- glm(defeitos ~ temperatura + densidade + taxa ,data = dados, 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),
data = dados)
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.4Conforme o resultado do método stepwise, levando em consideração o critério do AIC, ficamos com apenas temperatura e densidade como variáveis independentes no ajuste a ser feito.
fit <- glm(defeitos ~ temperatura + densidade, family = gaussian(link=identity), data = dados)
summary(fit)
Call:
glm(formula = defeitos ~ temperatura + densidade, family = gaussian(link = identity),
data = dados)
Deviance Residuals:
Min 1Q Median 3Q Max
-14.2233 -4.2245 -0.5605 3.8984 15.5638
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 46.238 52.066 0.888 0.3823
temperatura 18.050 7.965 2.266 0.0317 *
densidade -2.327 1.383 -1.683 0.1040
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 50.1779)
Null deviance: 10923.9 on 29 degrees of freedom
Residual deviance: 1354.8 on 27 degrees of freedom
AIC: 207.44
Number of Fisher Scoring iterations: 2O modelo ajustado, tem como equação:
\[defeitos_i = 46.24 + 18.05\times temperatura_i -2.33 \times densidade_i \]
Ao aumentarmos a temperatura do processo de produção, o número médio de defeitos nos produtos aumenta 18.05.
Se aumentarmos a quantidade da densidade do produto, diminuirá o número médio nos produtos em 2.327.
ajuste<-summary(fit)
knitr::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 |
- Realize uma análise de diagnóstico para o melhor modelo obtido em (a). Analise os resultados obtidos.
Análise Diagnóstico
Normalidade
rd <- boot::glm.diag(fit)$rd
teste1<-nortest::lillie.test(rd)
teste2<-shapiro.test(rd)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 0.4610346 | 0.5846342 |
Podemos observar pela tabela acima que os resíduos são normais. Também podemos verificar a normalidade através do gráfico qq-plot.
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample=rd))+
ggplot2::stat_qq(size=4, color="#e7ad52") + ggplot2::ylim(-4, 4) + ggplot2::xlim(-4, 4) +
ggplot2::stat_qq_line(colour = '#001132')+ ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
q1
É comprovado através dos testes e da análise gráfica que os resíduos são distribuidos normalmente.
Pontos aberrantes
indice=as.vector(1:length(defeitos))
rd1 <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=rd)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = -2, slope = 0, color='#001132')+ ggplot2::geom_abline(intercept = 2, slope = 0, color='#001132') + ggplot2::ylim(-4,4) + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
rd1
Visualmente, é observado que duas observações parecem estar fora dos limites estabelecidos, sendo assim dois possíveis pontos aberrante, no qual é a observação de número 4 e 28.
dados[which(rd>2 | rd < -2),] temperatura densidade taxa defeitos
4 2.36 26.3 222.1 13.4
28 2.82 22.3 253.2 60.8Pontos de Alavanca
h<- boot::glm.diag(fit)$h
h3 <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=h)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(fit))/length(defeitos)), slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
h3
Foi encontrado um ponto de alavanca, que é a observação de número 14. Sendo as observações deste ponto de alavanca:
dados[which(h>2*(length(coef(fit))/length(defeitos))),] temperatura densidade taxa defeitos
14 2.73 24.85 251.9 37.6Pontos de Influência
DC<-boot::glm.diag(fit)$cook
dc <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.05,length(coef(fit)))/2, slope = 0, color='#001132') +
ggplot2::ylim(0,0.5) + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
dc
Conforme observamos acima, há dois pontos de influência, que são as observações de número 1 e 4. Sendo os valores das observações:
dados[which(DC>qchisq(0.05,length(coef(fit)))/2),] temperatura densidade taxa defeitos
1 0.97 32.08 177.7 0.2
4 2.36 26.30 222.1 13.4