rm(list=ls(all=T))Lista de Exercícios
Questão 4
O conjunto de dados descrito no arquivo censo.dat, extraído do censo do IBGE de 2000, apresenta para cada unidade da federação o número médio de anos de estudo e a renda média mensal (em reais) do chefe ou chefes do domicílio. Suponha que a renda tenha distribuição gama.
Definindo o diretório
O diretório é importante, pois é onde se encontra os dados.
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/LISTA 1")Importando o banco de dados
dados <- read.table("censo.txt", header = F,
col.names = c("UF", "estudo", "renda"), dec = ".")
dados UF estudo renda
1 RR 5.7 685
2 AP 6.0 683
3 AC 4.5 526
4 RO 4.9 662
5 PA 4.7 536
6 AM 5.5 627
7 TO 4.5 520
8 PB 3.9 423
9 MA 3.6 343
10 RN 4.5 513
11 SE 4.3 462
12 PI 3.5 383
13 BA 4.1 460
14 PE 4.6 517
15 AL 3.7 454
16 CE 4.0 448
17 SP 6.8 1076
18 RJ 7.1 970
19 ES 5.7 722
20 MG 5.4 681
21 SC 6.3 814
22 RS 6.4 800
23 PR 6.0 782
24 MT 5.4 775
25 GO 5.5 689
26 MS 5.7 731
27 DF 8.2 1499summary(dados) UF estudo renda
Length:27 Min. :3.500 Min. : 343.0
Class :character 1st Qu.:4.400 1st Qu.: 487.5
Mode :character Median :5.400 Median : 662.0
Mean :5.204 Mean : 658.6
3rd Qu.:5.850 3rd Qu.: 753.0
Max. :8.200 Max. :1499.0 Histograma
O histograma tem a finalidade de ver como os dados estão distribuidos.
ggplot2::ggplot(dados, ggplot2::aes(x = dados$renda))+
ggplot2::geom_histogram(fill='white',
color = "blue",
breaks = hist(dados$renda, plot = F)$breaks) + ggplot2::xlab("Renda") + ggplot2::ylab("Frequência") + ggplot2::ggtitle("Histograma da Renda")Warning: Use of `dados$renda` is discouraged.
ℹ Use `renda` instead.
hist(dados$renda, xlab = "Renda", ylab = "Frequência", main = "Histograma da Renda")
Ajuste dos modelos com a escolha das funções de ligação
- Realize o ajuste da gama 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(renda~estudo, family = Gamma(link=identity))
fit2<-glm(renda~estudo, family = Gamma(link=log))
fit3<-glm(renda~estudo, family = Gamma(link=inverse))Análise de Desvio
Vamos verificar a adequação dos modelos. Se \(D(y;\mu)/\phi \leq \chi^{2}_{n-p;1-\alpha}\), 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)Verificando se todos os modelos foram aceitos:
| Modelo 1 | Modelo 2 | Modelo 3 | |
|---|---|---|---|
| Desvio/\(\phi\) | 23.5916959 | 24.625375 | 25.0564704 |
| \(\chi^2\) | 37.6524841 | 37.6524841 | 37.6524841 |
Pela tabela acima podemos verificar que todos os modelos são adequados.
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<-rendag1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va1,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(300,1500) + ggplot2::xlim(300,1200) + 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(300,1500) + ggplot2::xlim(300,1500) +
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(300,1500) + ggplot2::xlim(300,1800)+
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)
Escolha da função de variância
Gráfico de \(r_D^{\ast}\) 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)
Com base na função de variância e função de ligação, vemos que visualmente o melhor modelo ajustado é com a função de ligação inversa.
Então, escolhemos esse como o modelo a ser adotado.
- Selecione as variáveis. Realize uma análise residual para o melhor modelo obtido. O modelo é adequado? Por quê?
Como estamos no caso, de uma regressão simples, apenas com uma variável, não precisamos fazer a seleção de variáveis. Vamos ficar apenas com o modelo com função de ligação inversa, sendo este o modelo definido.
summary(fit3)
Call:
glm(formula = renda ~ estudo, family = Gamma(link = inverse))
Deviance Residuals:
Min 1Q Median 3Q Max
-0.25264 -0.05819 -0.01238 0.06872 0.21644
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.539e-03 1.243e-04 28.48 < 2e-16 ***
estudo -3.608e-04 1.971e-05 -18.31 5.37e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.01044499)
Null deviance: 3.03481 on 26 degrees of freedom
Residual deviance: 0.26171 on 25 degrees of freedom
AIC: 304.91
Number of Fisher Scoring iterations: 4ajuste<-summary(fit3)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 0.0035385 | 0.0001242 | 28.48024 | 0 |
| estudo | -0.0003608 | 0.0000197 | -18.31010 | 0 |
No qual o modelo de regressão ficou expresso como:
\[\frac{1}{\mu_i} = 0.00354 -3.6\times 10^{-4}\times estudo_i\]
Conforme temos o acréscimo de uma unidade na quantidade de anos estudados, vamos ter uma diminuição de \(3.6 \times 10^{-4}\) valor no inverso da renda média.
Análise de Diagnósticos
Teste de Normalidade
Testar se os resíduos são normais.
teste1<-nortest::lillie.test(rd3)
teste2<-shapiro.test(rd3)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 0.4438252 | 0.4438171 |
Podemos observar pela tabela acima que os resíduos seguem uma distribuição normal. Também podemos verificar a normalidade através do gráfico qq-plot.
q1 <- 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::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
O gráfico mais uma vez comprova a normalidade dos resíduos
- Realize uma análise de diagnóstico para o melhor modelo obtido em (a). Analise os resultados obtidos.
Diagnósticos
Aberrantes
indice=as.vector(1:length(renda))
rd <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=rd3)) + 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'))
rd
Conforme vemos a análise gráfica, temos 3 pontos aberrantes, que são os pontos 9, 24 e 27. Sendo as observações:
dados[which(rd3>2 | rd3< -2),] UF estudo renda
9 MA 3.6 343
24 MT 5.4 775
27 DF 8.2 1499Alavanca
h<- boot::glm.diag(fit3)$h
indice=as.vector(1:length(renda))
h3 <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=h)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(fit3))/length(renda)), 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
Conforme o gráfico vemos que tem um ponto de alavanca nos dados que é a observação de número 27.
dados[which(h>2*(length(coef(fit3))/length(renda))),] UF estudo renda
27 DF 8.2 1499Influentes
DC<-boot::glm.diag(fit3)$cook
dc <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.1,length(coef(fit3)))/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
Mais uma vez temos a presença de três pontos Influente, que são as observação de número: 9, 17 e 27. Que são as observações abaixo:
dados[which(DC>qchisq(0.1,length(coef(fit3)))/2),] UF estudo renda
9 MA 3.6 343
17 SP 6.8 1076
27 DF 8.2 1499Questão 5
- O conjunto de dados
DoctorVisitsna biblioteca AER do R, extraído do Australian Health Survey de 1977-78, com 5190 homens tanto jovens como velhos. A variável resposta évisits.
dados <- read.csv2("DoctorVisits.csv", sep = ",")Fazendo uma breve descritiva
str(dados)'data.frame': 5190 obs. of 13 variables:
$ X : int 1 2 3 4 5 6 7 8 9 10 ...
$ visits : int 1 1 1 1 1 1 1 1 1 1 ...
$ gender : chr "female" "female" "male" "male" ...
$ age : chr "0.19" "0.19" "0.19" "0.19" ...
$ income : chr "0.55" "0.45" "0.9" "0.15" ...
$ illness : int 1 1 3 1 2 5 4 3 2 1 ...
$ reduced : int 4 2 0 0 5 1 0 0 0 0 ...
$ health : int 1 1 0 0 1 9 2 6 5 0 ...
$ private : chr "yes" "yes" "no" "no" ...
$ freepoor : chr "no" "no" "no" "no" ...
$ freerepat: chr "no" "no" "no" "no" ...
$ nchronic : chr "no" "no" "no" "no" ...
$ lchronic : chr "no" "no" "no" "no" ...dados$age <- as.numeric(dados$age)*100
dados$income <- (as.numeric(dados$income)/12)*10000DescTools::Desc(dados, plotit = F, digits = 2)------------------------------------------------------------------------------
Describe dados (data.frame):
data frame: 5190 obs. of 13 variables
5190 complete cases (100.0%)
Nr ColName Class NAs Levels
1 X integer .
2 visits integer .
3 gender character .
4 age numeric .
5 income numeric .
6 illness integer .
7 reduced integer .
8 health integer .
9 private character .
10 freepoor character .
11 freerepat character .
12 nchronic character .
13 lchronic character .
------------------------------------------------------------------------------
1 - X (integer)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 = n 0 2'595.50 2'554.73
100.0% 0.0% 0.0% 2'636.27
.05 .10 .25 median .75 .90 .95
260.45 519.90 1'298.25 2'595.50 3'892.75 4'671.10 4'930.55
range sd vcoef mad IQR skew kurt
5'189.00 1'498.37 0.58 1'923.67 2'594.50 0.00 -1.20
lowest : 1, 2, 3, 4, 5
highest: 5'186, 5'187, 5'188, 5'189, 5'190
' 95%-CI (classic)
------------------------------------------------------------------------------
2 - visits (integer)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 10 4'141 0.30 0.28
100.0% 0.0% 79.8% 0.32
.05 .10 .25 median .75 .90 .95
0.00 0.00 0.00 0.00 0.00 1.00 2.00
range sd vcoef mad IQR skew kurt
9.00 0.80 2.65 0.00 0.00 4.74 31.29
value freq perc cumfreq cumperc
1 0 4'141 79.8% 4'141 79.8%
2 1 782 15.1% 4'923 94.9%
3 2 174 3.4% 5'097 98.2%
4 3 30 0.6% 5'127 98.8%
5 4 24 0.5% 5'151 99.2%
6 5 9 0.2% 5'160 99.4%
7 6 12 0.2% 5'172 99.7%
8 7 12 0.2% 5'184 99.9%
9 8 5 0.1% 5'189 100.0%
10 9 1 0.0% 5'190 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
3 - gender (character - dichotomous)
length n NAs unique
5'190 5'190 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
female 2'702 52.06% 50.70% 53.42%
male 2'488 47.94% 46.58% 49.30%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
4 - age (numeric)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 12 0 40.64 40.08
100.0% 0.0% 0.0% 41.20
.05 .10 .25 median .75 .90 .95
19.00 19.00 22.00 32.00 62.00 72.00 72.00
range sd vcoef mad IQR skew kurt
53.00 20.48 0.50 19.27 40.00 0.42 -1.50
value freq perc cumfreq cumperc
1 19 752 14.5% 752 14.5%
2 22 1'213 23.4% 1'965 37.9%
3 27 523 10.1% 2'488 47.9%
4 32 301 5.8% 2'789 53.7%
5 37 146 2.8% 2'935 56.6%
6 42 126 2.4% 3'061 59.0%
7 47 181 3.5% 3'242 62.5%
8 52 222 4.3% 3'464 66.7%
9 57 273 5.3% 3'737 72.0%
10 62 316 6.1% 4'053 78.1%
11 67 315 6.1% 4'368 84.2%
12 72 822 15.8% 5'190 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
5 - income (numeric)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 14 79 485.97 477.60
100.0% 0.0% 1.5% 494.33
.05 .10 .25 median .75 .90 .95
125.00 208.33 208.33 458.33 750.00 916.67 1'083.33
range sd vcoef mad IQR skew kurt
1'250.00 307.42 0.63 370.65 541.67 0.73 -0.18
lowest : 0.00 (79), 8.33 (35), 50.00 (80), 125.00 (249), 208.33 (1'195)
highest: 625.00 (441), 750.00 (589), 916.67 (361), 1'083.33 (162), 1'250.00 (215)
heap(?): remarkable frequency (23.0%) for the mode(s) (= 208.333333333333)
' 95%-CI (classic)
------------------------------------------------------------------------------
6 - illness (integer)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 6 1'554 1.43 1.39
100.0% 0.0% 29.9% 1.47
.05 .10 .25 median .75 .90 .95
0.00 0.00 0.00 1.00 2.00 3.00 4.00
range sd vcoef mad IQR skew kurt
5.00 1.38 0.97 1.48 2.00 0.94 0.16
value freq perc cumfreq cumperc
1 0 1'554 29.9% 1'554 29.9%
2 1 1'638 31.6% 3'192 61.5%
3 2 946 18.2% 4'138 79.7%
4 3 542 10.4% 4'680 90.2%
5 4 274 5.3% 4'954 95.5%
6 5 236 4.5% 5'190 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
7 - reduced (integer)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 15 4'454 0.86 0.78
100.0% 0.0% 85.8% 0.94
.05 .10 .25 median .75 .90 .95
0.00 0.00 0.00 0.00 0.00 2.00 7.00
range sd vcoef mad IQR skew kurt
14.00 2.89 3.35 0.00 0.00 3.83 13.82
lowest : 0 (4'454), 1 (177), 2 (108), 3 (74), 4 (45)
highest: 10 (12), 11 (2), 12 (6), 13 (5), 14 (188)
heap(?): remarkable frequency (85.8%) for the mode(s) (= 0)
' 95%-CI (classic)
------------------------------------------------------------------------------
8 - health (integer)
length n NAs unique 0s mean meanCI'
5'190 5'190 0 13 3'026 1.22 1.16
100.0% 0.0% 58.3% 1.28
.05 .10 .25 median .75 .90 .95
0.00 0.00 0.00 0.00 2.00 4.00 6.00
range sd vcoef mad IQR skew kurt
12.00 2.12 1.74 0.00 2.00 2.40 6.26
lowest : 0 (3'026), 1 (823), 2 (446), 3 (273), 4 (187)
highest: 8 (42), 9 (32), 10 (21), 11 (24), 12 (19)
heap(?): remarkable frequency (58.3%) for the mode(s) (= 0)
' 95%-CI (classic)
------------------------------------------------------------------------------
9 - private (character - dichotomous)
length n NAs unique
5'190 5'190 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
no 2'892 55.72% 54.37% 57.07%
yes 2'298 44.28% 42.93% 45.63%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
10 - freepoor (character - dichotomous)
length n NAs unique
5'190 5'190 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
no 4'968 95.72% 95.14% 96.24%
yes 222 4.28% 3.76% 4.86%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
11 - freerepat (character - dichotomous)
length n NAs unique
5'190 5'190 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
no 4'099 78.98% 77.85% 80.07%
yes 1'091 21.02% 19.93% 22.15%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
12 - nchronic (character - dichotomous)
length n NAs unique
5'190 5'190 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
no 3'098 59.69% 58.35% 61.02%
yes 2'092 40.31% 38.98% 41.65%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
13 - lchronic (character - dichotomous)
length n NAs unique
5'190 5'190 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
no 4'585 88.34% 87.44% 89.19%
yes 605 11.66% 10.81% 12.56%
' 95%-CI (Wilson)attach(dados)Histograma
O histograma tem a finalidade de ver como os dados estão distribuidos.
ggplot2::ggplot(dados, ggplot2::aes(x = visits))+
ggplot2::geom_histogram(fill='white',
color = "blue", breaks = hist(visits, plot = F)$breaks) + ggplot2::xlab("Número de Visitas") + ggplot2::ylab("Frequência") +
ggplot2::ggtitle("Histograma do número de visitas")
hist(visits, xlab = "Número de Visitas",
ylab = "Frequência", main = "Histograma do número de Visitas")
Como é observado no histograma e através da análise descritiva foi visto que os dados por serem inteiros, se aproximam e serão modelados através da distribuição Poisson.
Transformando em fatores
gender <- as.factor(gender)
private <- as.factor(private)
freepoor <- as.factor(freepoor)
freerepat <- as.factor(freerepat)
nchronic <- as.factor(nchronic)
lchronic <- as.factor(lchronic)- 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(visits ~ age + income + illness +
reduced + health + gender + private + freepoor +
freerepat + nchronic + lchronic,
family = poisson(link=log), data = dados)
fit2<-glm(visits ~ age + gender + income + illness +
reduced + health,
family = poisson(link=sqrt), data = dados)Vamos utilizar todas as variáveis nesse primeiro momento e depois na seleção de variáveis saberemos quais são as significantes que ficará no modelo final.
Análise de Desvio
Vamos verificar a adequação dos modelos. Se \(D(y;\mu)/\phi \leq \chi^{2}_{n-p;1-\alpha}\), 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)Verificando se todos os modelos foram aceitos:
| Modelo 1 | Modelo 2 | |
|---|---|---|
| Desvio/\(\phi\) | 4380.1331067 | 4214.1430688 |
| \(\chi^2\) | 5346.516891 | 5351.597692 |
Como é observado todos os dois 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)
y<-visitssummary(y) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0000 0.0000 0.3017 0.0000 9.0000 summary(va1) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.07314 0.15235 0.20196 0.30173 0.28978 4.35183 summary(va2) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.05323 0.10829 0.18354 0.30173 0.30664 3.04994 g1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va1,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,10) + ggplot2::xlim(0,6) + 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(dados, ggplot2::aes(x=va2,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,10) + ggplot2::xlim(0,6) +
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)
Escolha da função de variância
Gráfico de \(r_D^{\ast}\) versus valores ajustados.
rd1 <- boot::glm.diag(fit1)$rd
rd2 <- boot::glm.diag(fit2)$rdq1 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd1))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,6) + ggplot2::xlim(-4,6) +
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(dados, ggplot2::aes(sample = rd2))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,6) + ggplot2::xlim(-4,6) +
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 podemos observar a função de ligação Sqrt parece se adequar melhor aos nossos dados.
- Selecione as variáveis. Realize uma análise residual para o melhor modelo obtido. O modelo é adequado? Por quê?
Seleção das Variáveis
menor <- glm(visits ~ 1, family = poisson(link = sqrt), data = dados)
maior <- glm(visits ~ gender + age + income + illness + reduced +
health + private + freepoor + freerepat + nchronic +
lchronic ,data = dados, family = poisson(link = sqrt))
step(menor, scope = list(upper = maior), direction = "both")Start: AIC=7968.39
visits ~ 1
Df Deviance AIC
+ reduced 1 4557.4 6893.0
+ illness 1 5112.2 7447.8
+ health 1 5311.2 7646.8
+ lchronic 1 5468.7 7804.3
+ age 1 5470.2 7805.8
+ freerepat 1 5523.3 7858.8
+ gender 1 5566.2 7901.7
+ income 1 5570.5 7906.1
+ nchronic 1 5612.8 7948.3
+ freepoor 1 5615.6 7951.2
<none> 5634.8 7968.4
+ private 1 5634.1 7969.7
Step: AIC=6893.01
visits ~ reduced
Df Deviance AIC
+ illness 1 4306.6 6644.2
+ age 1 4449.7 6787.3
+ freerepat 1 4472.1 6809.7
+ health 1 4488.8 6826.4
+ gender 1 4497.7 6835.3
+ income 1 4514.1 6851.7
+ lchronic 1 4520.8 6858.4
+ nchronic 1 4525.3 6862.8
+ freepoor 1 4537.9 6875.4
<none> 4557.4 6893.0
+ private 1 4557.4 6894.9
- reduced 1 5634.8 7968.4
Step: AIC=6644.16
visits ~ reduced + illness
Df Deviance AIC
+ age 1 4251.8 6591.3
+ freerepat 1 4264.3 6603.8
+ gender 1 4269.4 6609.0
+ income 1 4287.2 6626.8
+ freepoor 1 4290.8 6630.4
+ health 1 4292.3 6631.9
+ lchronic 1 4297.2 6636.7
+ nchronic 1 4301.9 6641.5
<none> 4306.6 6644.2
+ private 1 4306.0 6645.5
- illness 1 4557.4 6893.0
- reduced 1 5112.2 7447.8
Step: AIC=6591.33
visits ~ reduced + illness + age
Df Deviance AIC
+ gender 1 4231.3 6572.9
+ health 1 4235.4 6576.9
+ freepoor 1 4244.3 6585.9
+ freerepat 1 4244.6 6586.2
+ lchronic 1 4245.0 6586.6
+ income 1 4246.0 6587.6
+ private 1 4249.6 6591.2
<none> 4251.8 6591.3
+ nchronic 1 4251.8 6593.3
- age 1 4306.6 6644.2
- illness 1 4449.7 6787.3
- reduced 1 5047.2 7384.7
Step: AIC=6572.85
visits ~ reduced + illness + age + gender
Df Deviance AIC
+ health 1 4215.7 6559.2
+ freepoor 1 4223.8 6567.4
+ lchronic 1 4224.2 6567.8
+ freerepat 1 4226.1 6569.6
+ income 1 4229.2 6572.8
<none> 4231.3 6572.9
+ private 1 4230.0 6573.5
+ nchronic 1 4231.2 6574.8
- gender 1 4251.8 6591.3
- age 1 4269.4 6609.0
- illness 1 4420.5 6760.1
- reduced 1 5029.5 7369.1
Step: AIC=6559.24
visits ~ reduced + illness + age + gender + health
Df Deviance AIC
+ freepoor 1 4207.3 6552.8
+ lchronic 1 4210.7 6556.2
+ freerepat 1 4211.8 6557.4
<none> 4215.7 6559.2
+ private 1 4214.0 6559.6
+ income 1 4214.1 6559.7
+ nchronic 1 4215.7 6561.2
- health 1 4231.3 6572.9
- gender 1 4235.4 6576.9
- age 1 4255.6 6597.2
- illness 1 4358.9 6700.5
- reduced 1 4911.1 7252.6
Step: AIC=6552.83
visits ~ reduced + illness + age + gender + health + freepoor
Df Deviance AIC
+ lchronic 1 4202.4 6549.9
+ freerepat 1 4203.5 6551.1
+ income 1 4203.5 6551.1
<none> 4207.3 6552.8
+ private 1 4206.7 6554.3
+ nchronic 1 4207.2 6554.8
- freepoor 1 4215.7 6559.2
- health 1 4223.8 6567.4
- gender 1 4226.9 6570.4
- age 1 4240.3 6583.9
- illness 1 4349.4 6693.0
- reduced 1 4902.7 7246.3
Step: AIC=6549.93
visits ~ reduced + illness + age + gender + health + freepoor +
lchronic
Df Deviance AIC
+ income 1 4198.6 6548.2
+ freerepat 1 4199.4 6549.0
<none> 4202.4 6549.9
+ nchronic 1 4201.0 6550.5
+ private 1 4201.8 6551.4
- lchronic 1 4207.3 6552.8
- freepoor 1 4210.7 6556.2
- health 1 4216.8 6562.4
- gender 1 4222.2 6567.8
- age 1 4233.6 6579.2
- illness 1 4334.8 6680.3
- reduced 1 4872.9 7218.5
Step: AIC=6548.21
visits ~ reduced + illness + age + gender + health + freepoor +
lchronic + income
Df Deviance AIC
<none> 4198.6 6548.2
+ private 1 4196.9 6548.4
+ nchronic 1 4197.1 6548.6
+ freerepat 1 4197.2 6548.8
- income 1 4202.4 6549.9
- lchronic 1 4203.5 6551.1
- freepoor 1 4209.2 6556.7
- health 1 4212.5 6560.0
- gender 1 4214.0 6561.6
- age 1 4222.6 6570.2
- illness 1 4328.4 6676.0
- reduced 1 4869.0 7216.6
Call: glm(formula = visits ~ reduced + illness + age + gender + health +
freepoor + lchronic + income, family = poisson(link = sqrt),
data = dados)
Coefficients:
(Intercept) reduced illness age gendermale health
3.044e-01 6.091e-02 6.480e-02 1.822e-03 -5.843e-02 1.292e-02
freepooryes lchronicyes income
-1.213e-01 4.885e-02 -4.731e-05
Degrees of Freedom: 5189 Total (i.e. Null); 5181 Residual
Null Deviance: 5635
Residual Deviance: 4199 AIC: 6548No final ficamos com essas 10 variáveis independentes no modelo ajustado:
O ajuste final expresso incorporou todas as variáveis.
fit <- glm(visits ~ reduced + illness + freerepat + gender +
private + health + freepoor + income + lchronic + nchronic,
family = poisson(link = sqrt), data = dados)
summary(fit)
Call:
glm(formula = visits ~ reduced + illness + freerepat + gender +
private + health + freepoor + income + lchronic + nchronic,
family = poisson(link = sqrt), data = dados)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4389 -0.6819 -0.5109 -0.3548 5.3815
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.235e-01 2.177e-02 14.862 < 2e-16 ***
reduced 6.155e-02 2.558e-03 24.063 < 2e-16 ***
illness 6.361e-02 5.839e-03 10.894 < 2e-16 ***
freerepatyes 9.513e-02 2.232e-02 4.262 2.03e-05 ***
gendermale -5.217e-02 1.507e-02 -3.463 0.000535 ***
privateyes 5.368e-02 1.686e-02 3.183 0.001457 **
health 1.203e-02 3.620e-03 3.323 0.000891 ***
freepooryes -1.117e-01 3.670e-02 -3.043 0.002345 **
income -5.459e-05 2.614e-05 -2.089 0.036748 *
lchronicyes 5.835e-02 2.477e-02 2.356 0.018496 *
nchronicyes 2.924e-02 1.619e-02 1.806 0.070877 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 5634.8 on 5189 degrees of freedom
Residual deviance: 4195.4 on 5179 degrees of freedom
AIC: 6548.9
Number of Fisher Scoring iterations: 5Podemos observar que o modelo final ficou como sendo expresso:
ajuste<-summary(fit)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 0.3235290 | 0.0217686 | 14.862210 | 0.0000000 |
| reduced | 0.0615500 | 0.0025578 | 24.063488 | 0.0000000 |
| illness | 0.0636127 | 0.0058395 | 10.893549 | 0.0000000 |
| freerepatyes | 0.0951267 | 0.0223203 | 4.261900 | 0.0000203 |
| gendermale | -0.0521682 | 0.0150664 | -3.462550 | 0.0005351 |
| privateyes | 0.0536805 | 0.0168640 | 3.183149 | 0.0014568 |
| health | 0.0120280 | 0.0036198 | 3.322823 | 0.0008911 |
| freepooryes | -0.1116661 | 0.0366997 | -3.042703 | 0.0023446 |
| income | -0.0000546 | 0.0000261 | -2.088558 | 0.0367476 |
| lchronicyes | 0.0583468 | 0.0247700 | 2.355541 | 0.0184958 |
| nchronicyes | 0.0292444 | 0.0161905 | 1.806264 | 0.0708772 |
As variáveis categoricas foram incorporadas como variáveis dummy, sendo significantes para o modelo final.
Análise de Diagnósticos
Teste de Normalidade
Testar se os resíduos são normais.
rd <- boot::glm.diag(fit)$rd
teste1<-nortest::lillie.test(rd)| Teste 1 | |
|---|---|
| p-valor | 0 |
Podemos observar pela tabela acima que os resíduos não são normais. Também podemos verificar a não 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, 6) + ggplot2::xlim(-4, 6) +
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
Conforme a análise gráfica acima é observado que os resíduos não seguem distribuição normal.
- Realize uma análise de diagnóstico para o melhor modelo obtido em (a). Analise os resultados obtidos.
Diagnósticos
Aberrantes
indice=as.vector(1:length(visits))
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(-6,6) + 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
De acordo com a análise gráfica temos várias observações como outlier. Prosseguiremos com a análise de diagnóstico. Mais precisamente, temos o número de 171.
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(visits)), 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
Como pontos de alavanca, temos um quantitativo de 510 pontos de alavanca.
Influentes
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.1,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')) + ggplot2::ylim(0,qchisq(0.1,length(coef(fit)))/2 + 1)
dc
Foi diagnosticado que não existe nenhum ponto de influência.
Questão 6
- No arquivo
vendas.datsão descritas informações a respeito das vendas no ano anterior de um tipo de telhado de madeira em 26 filiais de uma rede de lojas de construção. As variáveis estão colocadas na seguinte ordem:telhados, total de telhados vendidos (em mil metros quadrados),gastos, gastos pela loja com promoções do produto (em mil USD),clientes, número de clientes cadastrados na loja (em milhares),marcas, número de marcas concorrentes do produto epotencial, potencial da loja (quanto maior o valor maior o potencial). A variável resposta ételhado.
rm(list = ls(all=T))
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/LISTA 1")dados <- read.table("vendas.txt", dec = ".", col.names = c("telhado", "gastos", "clientes", "marcas", "potencial"))
stargazer::stargazer(dados, summary = F, type = "text")
===========================================
telhado gastos clientes marcas potencial
-------------------------------------------
1 79.300 5.500 31 10 8
2 200.100 2.500 55 8 6
3 163.200 8 67 12 9
4 200.100 3 50 7 16
5 146 3 38 8 15
6 177.700 2.900 71 12 17
7 30.900 8 30 12 8
8 291.900 9 56 5 10
9 160 4 42 8 4
10 339.400 6.500 73 5 16
11 159.600 5.500 60 11 7
12 86.300 5 44 12 12
13 237.500 6 50 6 6
14 107.200 5 39 10 4
15 155 3.500 55 10 4
16 291.400 8 70 6 14
17 100.200 6 40 11 6
18 135.800 4 50 11 8
19 223.300 7.500 62 9 13
20 195 7 59 9 11
21 73.400 6.700 53 13 5
22 47.700 6.100 38 13 10
23 140.700 3.600 43 9 17
24 93.500 4.200 26 8 3
25 259 4.500 75 8 19
26 331.200 5.600 71 4 9
-------------------------------------------stargazer::stargazer(dados, summary = T, type = "text")
============================================
Statistic N Mean St. Dev. Min Max
--------------------------------------------
telhado 26 170.208 84.549 30.900 339.400
gastos 26 5.408 1.832 2.500 9.000
clientes 26 51.846 14.215 26 75
marcas 26 9.115 2.582 4 13
potencial 26 9.885 4.727 3 19
--------------------------------------------Histograma
O histograma tem a finalidade de ver como os dados estão distribuidos.
attach(dados)
ggplot2::ggplot(dados, ggplot2::aes(x = telhado))+
ggplot2::geom_histogram(fill='white', color = "blue", breaks = hist(telhado, plot = F)$breaks) + ggplot2::xlab("Venda de Telhados") + ggplot2::ylab("Frequência") + ggplot2::ggtitle("Histograma da venda de telhados")
hist(telhado, xlab = "Venda de telhados", ylab = "Frequência", main = "Histograma da venda de telhados")
Podemos ver que os dados se adequam a uma distribuição normal.
- 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.
fit1<-glm(telhado~gastos + clientes + marcas + potencial, family = gaussian(link=identity), data = dados)
fit2<-glm(telhado~gastos + clientes + marcas + potencial, family = gaussian(link=log), data = dados)
fit3<-glm(telhado~gastos + clientes + marcas + potencial, family = gaussian(link=inverse), data = dados)Análise de Desvio
Vamos verificar a adequação dos modelos. Se \(D(y;\mu)/\phi \leq \chi^{2}_{n-p;1-\alpha}\), 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)Verificando se todos os modelos foram aceitos:
| Modelo 1 | Modelo 2 | Modelo 3 | |
|---|---|---|---|
| Desvio/\(\phi\) | 21 | 20.9996321 | 21.0007211 |
| \(\chi^2\) | 32.6705733 | 32.6705733 | 32.6705733 |
Pela tabela acima podemos verificar que todos os modelos são adequados.
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<-telhadog1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va1,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,400) + ggplot2::xlim(0,400) + 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,400) + ggplot2::xlim(0,400) +
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,400) + ggplot2::xlim(0,400)+
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)
Visualmente a função de ligação identidade se comporta melhor em nosso caso.
Escolha da função de variância
Gráfico de \(r_D^{\ast}\) 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)
Com base na função de variância, vemos que visualmente o melhor modelo ajustado é com a função de ligação identidade.
Então, escolhemos esse como o modelo a ser adotado.
- Selecione as variáveis. Realize uma análise residual para o melhor modelo obtido. O modelo é adequado? Por quê?
Seleção das Variáveis
menor <- glm(telhado~1, family = gaussian(link = identity), data = dados)
maior <- glm(telhado~gastos + clientes + marcas + potencial ,data = dados, family = gaussian(link = identity))
step(menor, scope = list(upper = maior), direction = "both")Start: AIC=307.51
telhado ~ 1
Df Deviance AIC
+ marcas 1 54712 278.73
+ clientes 1 69190 284.83
+ potencial 1 149071 304.79
<none> 178714 307.51
+ gastos 1 174204 308.84
Step: AIC=278.73
telhado ~ marcas
Df Deviance AIC
+ clientes 1 2210 197.30
+ potencial 1 44073 275.11
<none> 54712 278.73
+ gastos 1 51825 279.32
- marcas 1 178714 307.51
Step: AIC=197.3
telhado ~ marcas + clientes
Df Deviance AIC
+ gastos 1 1982 196.46
<none> 2210 197.30
+ potencial 1 2195 199.12
- clientes 1 54712 278.73
- marcas 1 69190 284.83
Step: AIC=196.46
telhado ~ marcas + clientes + gastos
Df Deviance AIC
<none> 1982 196.46
- gastos 1 2210 197.30
+ potencial 1 1937 197.87
- clientes 1 51825 279.32
- marcas 1 69086 286.80
Call: glm(formula = telhado ~ marcas + clientes + gastos, family = gaussian(link = identity),
data = dados)
Coefficients:
(Intercept) marcas clientes gastos
179.844 -21.217 3.369 1.677
Degrees of Freedom: 25 Total (i.e. Null); 22 Residual
Null Deviance: 178700
Residual Deviance: 1982 AIC: 196.5Podemos observar que o modelo final ficou como sendo expresso:
fit <- glm(telhado~marcas+clientes+gastos,family = gaussian(link = identity), data = dados)
summary(fit)
Call:
glm(formula = telhado ~ marcas + clientes + gastos, family = gaussian(link = identity),
data = dados)
Deviance Residuals:
Min 1Q Median 3Q Max
-20.4450 -4.6552 0.6802 6.1668 14.3571
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 179.8443 12.6213 14.249 1.37e-12 ***
marcas -21.2165 0.7773 -27.295 < 2e-16 ***
clientes 3.3694 0.1432 23.524 < 2e-16 ***
gastos 1.6772 1.0521 1.594 0.125
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 90.0697)
Null deviance: 178714.2 on 25 degrees of freedom
Residual deviance: 1981.5 on 22 degrees of freedom
AIC: 196.46
Number of Fisher Scoring iterations: 2De acordo com o resultado, foi obtido que a variável gastos não é tão significativa, pois o seu \(\beta\) existe evidências de que não é significativo.
Todavia, obtivemos que o modelo ajustado, tem como equação:
\[telhado_i = 179.84 -21.22\times marcas_i + 3.37 \times clientes_i + 1.68 \times gastos_i\]
Ao aumentarmos os gastos pela loja com promoções do produto (em mil USD) e o número de clientes cadastrados na loja (em milhares), o total de telhados vendidos aumenta em mil 1.677 metros quadrados e 3.369 mil metros quadrados, respectivamente.
Se aumentarmos o número de marcas concorrentes do produto, diminuirá o total de telhados vendidos em 21.217 mil metros quadrados.
ajuste<-summary(fit)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 179.844285 | 12.6213290 | 14.249235 | 0.0000000 |
| marcas | -21.216510 | 0.7772988 | -27.295179 | 0.0000000 |
| clientes | 3.369392 | 0.1432307 | 23.524241 | 0.0000000 |
| gastos | 1.677243 | 1.0520933 | 1.594196 | 0.1251589 |
Análise de Diagnósticos
Teste de Normalidade
Testar se os resíduos são normais.
rd <- boot::glm.diag(fit)$rd
teste1<-nortest::lillie.test(rd)
teste2<-shapiro.test(rd)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 0.9216333 | 0.8598839 |
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.
- Realize uma análise de diagnóstico para o melhor modelo obtido em (a). Analise os resultados obtidos.
Diagnósticos
Aberrantes
indice=as.vector(1:length(telhado))
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
De acordo com a análise gráfica temos apenas uma observação apontada como aberrante.
Que é a observação
dados[which(rd>2 | rd< - 2),] telhado gastos clientes marcas potencial
21 73.4 6.7 53 13 5Alavanca
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(telhado)), 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
Podemos ver um único ponto de alavanca, que é a observação de número 6. E os valores são os observados abaixo:
dados[which(h>2*(length(coef(fit))/length(telhado))),] telhado gastos clientes marcas potencial
6 177.7 2.9 71 12 17Influentes
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.1,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')) + ggplot2::ylim(-0.25,0.5)
dc
De acordo com a análise gráfica não temos nenhum ponto de influência.
Questão 7
- Em um experimento, rodas de turbinas foram usados por um número de horas, e o número de fissuras foram contadas. No arquivo
turbinas.datsão descritas as variáveis:horas, número de horas que a turbina foi usada,turbina, número de turbinas usadas para uma dada hora efissuras, número de fissuras nas turbina (variável resposta).
dados <- read.table("turbinas.txt", sep = ",", col.names = c("horas", "turbinas", "fissuras"))
stargazer::stargazer(dados, summary = F, type = "text")
==========================
horas turbinas fissuras
--------------------------
1 400 39 0
2 1,000 53 4
3 1,400 33 2
4 1,800 73 7
5 2,200 30 5
6 2,600 39 9
7 3,000 42 9
8 3,400 13 6
9 3,800 34 22
10 4,200 40 21
11 4,600 36 21
--------------------------summary(dados) horas turbinas fissuras
Min. : 400 Min. :13.00 Min. : 0.000
1st Qu.:1600 1st Qu.:33.50 1st Qu.: 4.500
Median :2600 Median :39.00 Median : 7.000
Mean :2582 Mean :39.27 Mean : 9.636
3rd Qu.:3600 3rd Qu.:41.00 3rd Qu.:15.000
Max. :4600 Max. :73.00 Max. :22.000 Histograma
O histograma tem a finalidade de ver como os dados estão distribuidos.
attach(dados)
ggplot2::ggplot(dados, ggplot2::aes(x = fissuras))+
ggplot2::geom_histogram(fill='white', color = "blue", breaks = hist(fissuras, plot = F)$breaks) + ggplot2::xlab("Quantidade de Fissuras") + ggplot2::ylab("Frequência") + ggplot2::ggtitle("Histograma da Quantidade de Fissuras")
hist(fissuras, xlab = "Quantidade de fissuras", ylab = "Frequência", main = "Histograma da quantidade de Fissuras")
Vamos modelar através de uma distribuição Poisson.
- 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(fissuras ~ horas + turbinas, family = poisson(link=log), data = dados)
fit2<-glm(fissuras~horas + turbinas, family = poisson(link=sqrt), data = dados)Análise de Desvio
Vamos verificar a adequação dos modelos. Se \(D(y;\mu)/\phi \leq \chi^{2}_{n-p;1-\alpha}\), 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)Verificando se todos os modelos foram aceitos:
| Modelo 1 | Modelo 2 | |
|---|---|---|
| Desvio/\(\phi\) | 8.984185 | 5.8429485 |
| \(\chi^2\) | 15.5073131 | 15.5073131 |
Como é observado todos os dois 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)
y<-fissurasg1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va1,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,30) + ggplot2::xlim(0,30) + 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(dados, ggplot2::aes(x=va2,y=y))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,30) + ggplot2::xlim(0,30) +
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)
Conforme podemos observar a função de ligação Raiz quadrada parece se adequar melhor aos nossos dados.
Escolha da função de variância
Gráfico de \(r_D^{\ast}\) versus valores ajustados.
rd1 <- boot::glm.diag(fit1)$rd
rd2 <- boot::glm.diag(fit2)$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("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(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("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)
Como é observado a escolha da função de variância o melhor modelo é com a função de ligação sqrt, então será este que adotaremos.
- Selecione as variáveis. Realize uma análise residual para o melhor modelo obtido. O modelo é adequado? Por quê?
Seleção das Variáveis
menor <- glm(fissuras ~ 1, family = poisson(link = sqrt), data = dados)
maior <- glm(fissuras ~ horas + turbinas ,data = dados, family = poisson(link = sqrt))
step(menor, scope = list(upper = maior), direction = "both")Start: AIC=110.9
fissuras ~ 1
Df Deviance AIC
+ horas 1 12.407 56.057
<none> 69.252 110.901
+ turbinas 1 68.953 112.603
Step: AIC=56.06
fissuras ~ horas
Df Deviance AIC
+ turbinas 1 5.843 51.492
<none> 12.407 56.057
- horas 1 69.252 110.901
Step: AIC=51.49
fissuras ~ horas + turbinas
Df Deviance AIC
<none> 5.843 51.492
- turbinas 1 12.407 56.057
- horas 1 68.953 112.603
Call: glm(formula = fissuras ~ horas + turbinas, family = poisson(link = sqrt),
data = dados)
Coefficients:
(Intercept) horas turbinas
-1.102107 0.001057 0.030477
Degrees of Freedom: 10 Total (i.e. Null); 8 Residual
Null Deviance: 69.25
Residual Deviance: 5.843 AIC: 51.49fit <- glm(fissuras~horas+turbinas, family = poisson(link = sqrt), data = dados)
summary(fit)
Call:
glm(formula = fissuras ~ horas + turbinas, family = poisson(link = sqrt),
data = dados)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.85627 -0.70420 0.04871 0.26391 1.51933
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.1021072 0.6500465 -1.695 0.08999 .
horas 0.0010574 0.0001239 8.532 < 2e-16 ***
turbinas 0.0304773 0.0113759 2.679 0.00738 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 69.2517 on 10 degrees of freedom
Residual deviance: 5.8429 on 8 degrees of freedom
AIC: 51.492
Number of Fisher Scoring iterations: 6O ajuste final expresso incorporou todas as variáveis, e ficou como é descrito abaixo:
\[\sqrt{\mu_i} = -1.1021 + 0.0011 \times horas_i + 0.0305 \times turbinas_i \]
Ao aumentarmos o número de horas que a turbina foi usada e o número de turbinas usadas para uma dada hora, a raiz quadrada do número esperado de fissuras nas turbina aumenta em 0.001 e 0.03 respectivamente.
ajuste<-summary(fit)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -1.1021072 | 0.6500465 | -1.695428 | 0.0899942 |
| horas | 0.0010574 | 0.0001239 | 8.532344 | 0.0000000 |
| turbinas | 0.0304773 | 0.0113759 | 2.679124 | 0.0073815 |
Análise de Diagnósticos
Teste de Normalidade
Testar se os resíduos são normais.
rd <- boot::glm.diag(fit)$rd
teste1<-nortest::lillie.test(rd)
teste2<-shapiro.test(rd)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 0.3943936 | 0.3964191 |
Podemos observar pela tabela acima que os resíduos estão distribuidos normalmente. 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 normais.
- Realize uma análise de diagnóstico para o melhor modelo obtido em (a). Analise os resultados obtidos.
Diagnósticos
Aberrantes
indice=as.vector(1:length(fissuras))
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
De acordo com a análise gráfica não temos nenhuma informação como outlier. Prosseguiremos com a análise de diagnóstico.
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(fissuras)), 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 diagnosticado um único ponto de alavanca, sendo ele a observação de número 4. Sendo os valores, os encontrados abaixo:
dados[which(h>2*(length(coef(fit))/length(fissuras))),] horas turbinas fissuras
4 1800 73 7Influentes
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.1,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
Foi observado um único ponto de influência que foi a observação de número 4.
Sendo a observação:
dados[which(DC>qchisq(0.1,length(coef(fit)))/2),] horas turbinas fissuras
4 1800 73 7