Prova 2 - Regressão II

Autor

Paulo Manoel da Silva Junior

Questão 1

  • O conjunto de dados descrito no arquivo heartdis.txt apresenta as variáveis caso, número do caso (desconsidere esta variável no modelo proposto) x1, pressão sistólica do sangue, x2, uma medida de colesterol, x3, variável dummy = 1 se há histórico na família de doenças cardíacas, x4, uma medida de obesidade, x5, idade e HeartDisease, se o paciente tem doençaa cardíaca (variável resposta).

Definindo o diretório

setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/PROVA 2")

Carregando o banco de dados

banco <- read.table("heartdis.txt", sep = ",", dec = ".", col.names = c("caso", "pressão_sistólica", "colesterol","dummy", "obesidade", "idade", "HeartDisease"), header = T)

Verificando a categoria das variáveis no banco e removendo a primeira coluna

banco <- banco[,-1]
glimpse(banco)
Rows: 462
Columns: 6
$ pressão_sistólica <int> 160, 144, 118, 170, 134, 132, 142, 114, 114, 132, 20…
$ colesterol        <dbl> 5.73, 4.41, 3.48, 6.41, 3.50, 6.47, 3.38, 4.59, 3.83…
$ dummy             <int> 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1…
$ obesidade         <dbl> 25.30, 28.87, 29.14, 31.99, 25.99, 30.77, 20.81, 23.…
$ idade             <int> 52, 63, 46, 58, 49, 45, 38, 58, 29, 53, 60, 40, 17, …
$ HeartDisease      <int> 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1…
banco$dummy <- as.factor(banco$dummy)
banco$HeartDisease <- as.factor(banco$HeartDisease)
glimpse(banco)
Rows: 462
Columns: 6
$ pressão_sistólica <int> 160, 144, 118, 170, 134, 132, 142, 114, 114, 132, 20…
$ colesterol        <dbl> 5.73, 4.41, 3.48, 6.41, 3.50, 6.47, 3.38, 4.59, 3.83…
$ dummy             <fct> 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1…
$ obesidade         <dbl> 25.30, 28.87, 29.14, 31.99, 25.99, 30.77, 20.81, 23.…
$ idade             <int> 52, 63, 46, 58, 49, 45, 38, 58, 29, 53, 60, 40, 17, …
$ HeartDisease      <fct> 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1…

Estatística Descritiva

Visualizando as primeiras linhas do banco

banco %>%
  slice_head(n=10) %>%
  gt::gt()
pressão_sistólica colesterol dummy obesidade idade HeartDisease
160 5.73 1 25.30 52 1
144 4.41 0 28.87 63 1
118 3.48 1 29.14 46 0
170 6.41 1 31.99 58 1
134 3.50 1 25.99 49 1
132 6.47 1 30.77 45 0
142 3.38 0 20.81 38 0
114 4.59 1 23.11 58 1
114 3.83 1 24.86 29 0
132 5.80 1 30.11 53 1
 
skimr::skim(banco)
Data summary
Name banco
Number of rows 462
Number of columns 6
_______________________
Column type frequency:
factor 2
numeric 4
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
dummy 0 1 FALSE 2 0: 270, 1: 192
HeartDisease 0 1 FALSE 2 0: 302, 1: 160

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
pressão_sistólica 0 1 138.33 20.50 101.00 124.00 134.00 148.00 218.00 ▅▇▃▁▁
colesterol 0 1 4.74 2.07 0.98 3.28 4.34 5.79 15.33 ▆▇▂▁▁
obesidade 0 1 26.04 4.21 14.70 22.98 25.80 28.50 46.58 ▂▇▃▁▁
idade 0 1 42.82 14.61 15.00 31.00 45.00 55.00 64.00 ▅▅▅▆▇

Letra A

  • Realize o ajuste da regressão logística e selecione as variáveis. O modelo é adequado? Justifique sua escolha.
menor <- glm(HeartDisease~1, data = banco, family = binomial(link="logit"))
maior <- glm(HeartDisease~pressão_sistólica+colesterol+obesidade+idade+dummy, data = banco, family = binomial(link = "logit"))
step(menor, scope = list(upper = maior), direction = "both")
Start:  AIC=598.11
HeartDisease ~ 1

                    Df Deviance    AIC
+ idade              1   525.56 529.56
+ dummy              1   561.89 565.89
+ colesterol         1   564.28 568.28
+ pressão_sistólica  1   579.32 583.32
+ obesidade          1   591.53 595.53
<none>                   596.11 598.11

Step:  AIC=529.56
HeartDisease ~ idade

                    Df Deviance    AIC
+ dummy              1   506.66 512.66
+ colesterol         1   512.48 518.48
<none>                   525.56 529.56
+ pressão_sistólica  1   524.44 530.44
+ obesidade          1   525.55 531.55
- idade              1   596.11 598.11

Step:  AIC=512.66
HeartDisease ~ idade + dummy

                    Df Deviance    AIC
+ colesterol         1   496.18 504.18
<none>                   506.66 512.66
+ pressão_sistólica  1   505.24 513.24
+ obesidade          1   506.63 514.63
- dummy              1   525.56 529.56
- idade              1   561.89 565.89

Step:  AIC=504.18
HeartDisease ~ idade + dummy + colesterol

                    Df Deviance    AIC
<none>                   496.18 504.18
+ obesidade          1   494.99 504.99
+ pressão_sistólica  1   495.17 505.17
- colesterol         1   506.66 512.66
- dummy              1   512.48 518.48
- idade              1   537.91 543.91

Call:  glm(formula = HeartDisease ~ idade + dummy + colesterol, family = binomial(link = "logit"), 
    data = banco)

Coefficients:
(Intercept)        idade       dummy1   colesterol  
   -4.35183      0.05476      0.88199      0.16980  

Degrees of Freedom: 461 Total (i.e. Null);  458 Residual
Null Deviance:      596.1 
Residual Deviance: 496.2    AIC: 504.2

Resultado: Após o procedimento da seleção de variáveis é observado que o modelo final ficou apenas com as seguintes variáveis: idade, sendo a idade do paciente, colesterol, conforme descrito acima que é uma medida de colesterol e dummy, que é igual a 1, se existe na família do paciente histórico de doenças cardíacas.

Ajustando o modelo e verificando se é adequado

fit <- glm(HeartDisease~idade+dummy+colesterol, data = banco, family = binomial(link="logit"))
phi11 <- summary(fit)$dispersion
desvio11 <- summary(fit)$deviance/phi11
q.quadr11 <- qchisq(0.95, desvio11)
Modelo
Desvio/\phi 496.1802755
\chi^2 549.107924

Temos do resultado acima que: O modelo em análise foi aceito.

ajuste<-summary(fit)
knitr::kable(ajuste$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.3518333 0.4912572 -8.858564 0.0000000
idade 0.0547551 0.0090766 6.032564 0.0000000
dummy1 0.8819917 0.2194687 4.018758 0.0000585
colesterol 0.1697957 0.0534461 3.176950 0.0014883

A equação é dada por:

log \left (\frac{\pi_1}{1-\pi_1} \right) = -4.352 + 0.055 \hspace{0.1cm}. \hspace{0.1cm}idade_i + 0.882\hspace{0.1cm} . \hspace{0.1cm}histórico_i + 0.1698\hspace{0.1cm} . \hspace{0.1cm}colesterol_i Resposta final: Podemos observar que o modelo ajustado é adequado, e que todas as variáveis independentes foram importantes para o modelo, pois no teste em que a hipótese nula é de que \beta_i = 0, essa hipótese foi rejeitada para todas as variáveis independentes do modelo ajustado.

Letra B

  • Faça a curva ROC do modelo. O que você pode concluir sobre o ajuste do modelo?
Epi::ROC(form = HeartDisease~idade+dummy+colesterol, plot = "ROC", MX = TRUE, data = banco)

Resposta: A curva ROC fornece que o modelo ficou bem ajustado, pois a suavização da curva fica evidente através da análise gráfica.

Letra C

  • Construa um envelope para os resíduos. Há algum ponto que não pertence ao envelope? Se sim, qual(is)?
par(mfrow=c(1,1))
X <- model.matrix(fit)
n <- nrow(X)
p <- ncol(X)
w <- fit$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
td <- resid(fit,type="deviance")/sqrt(1-h)
e <- matrix(0,n,100)
#
for(i in 1:100){
dif <- runif(n) - fitted(fit)
dif[dif >= 0 ] <- 0
dif[dif<0] <- 1
nresp <- dif
fit <- glm(nresp ~ X, family=binomial)
w <- fit$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
#
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:n){
  eo <- sort(e[i,])
e1[i] <- (eo[2]+eo[3])/2
e2[i] <- (eo[97]+eo[98])/2}
#
med <- apply(e,1,mean)
faixa <- range(td,e1,e2)
par(pty="s")
qqnorm(td,xlab="Percentil da N(0,1)",
ylab="Componente do Desvio", ylim=faixa, pch=16)
#
par(new=T)
#
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2)

Resposta: Não, aparentemente não existe nenhum ponto fora do envelope.

Letra D

  • Construa um intervalo de confiança de 90% para os parâmetros do modelo.

Na construção do envelope o objeto que tinha sido modelado acima foi renomeado, por isso precisamos ajustar o mesmo modelo só que com um nome diferente.

fit2 <- glm(HeartDisease~idade+dummy+colesterol, data = banco, family = binomial(link="logit"))
knitr::kable(confint(fit2,level = 0.90))
Waiting for profiling to be done...
5 % 95 %
(Intercept) -5.1911111 -3.5725014
idade 0.0401443 0.0700499
dummy1 0.5221366 1.2447973
colesterol 0.0830075 0.2593495

Letra E

  • Interprete o coeficiente \beta_5 da idade. Mantendo-se as outras variáveis constantes, o acréscimo de um ano na idade do paciente aumenta (ou diminui) em quanto a chance do paciente desenvolver uma doença cardíaca?

Resposta: Como foi analisado no ajuste acima, com o aumento de uma unidade da idade a chance do paciente desdenvolver uma doença cardíaca aumenta em 0.0548.

Questão 2

  • Considere o banco de dados Prestige do pacote carData do R que fornece 102 observações com seis variáveis das quais iremos utilizar apenas as variáveis: prestige (variável resposta) score de prestígio de Pineo-Porter para a ocupação, de uma pesquisa social feita nos meados dos anos 60, income, renda média, em dólares em 1971 e education, média, em anos, de estudo para a determinada educação.

Verificando a categoria das variáveis no banco

glimpse(carData::Prestige)
Rows: 102
Columns: 6
$ education <dbl> 13.11, 12.26, 12.77, 11.42, 14.62, 15.64, 15.09, 15.44, 14.5…
$ income    <int> 12351, 25879, 9271, 8865, 8403, 11030, 8258, 14163, 11377, 1…
$ women     <dbl> 11.16, 4.02, 15.70, 9.11, 11.68, 5.13, 25.65, 2.69, 1.03, 0.…
$ prestige  <dbl> 68.8, 69.1, 63.4, 56.8, 73.5, 77.6, 72.6, 78.1, 73.1, 68.8, …
$ census    <int> 1113, 1130, 1171, 1175, 2111, 2113, 2133, 2141, 2143, 2153, …
$ type      <fct> prof, prof, prof, prof, prof, prof, prof, prof, prof, prof, …

Estatística Descritiva

Visualizando as primeiras linhas do banco

carData::Prestige %>%
  slice_head(n=10) %>%
  gt::gt()
education income women prestige census type
13.11 12351 11.16 68.8 1113 prof
12.26 25879 4.02 69.1 1130 prof
12.77 9271 15.70 63.4 1171 prof
11.42 8865 9.11 56.8 1175 prof
14.62 8403 11.68 73.5 2111 prof
15.64 11030 5.13 77.6 2113 prof
15.09 8258 25.65 72.6 2133 prof
15.44 14163 2.69 78.1 2141 prof
14.52 11377 1.03 73.1 2143 prof
14.64 11023 0.94 68.8 2153 prof
 
skimr::skim(carData::Prestige)
Data summary
Name carData::Prestige
Number of rows 102
Number of columns 6
_______________________
Column type frequency:
factor 1
numeric 5
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
type 4 0.96 FALSE 3 bc: 44, pro: 31, wc: 23

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
education 0 1 10.74 2.73 6.38 8.45 10.54 12.65 15.97 ▆▇▆▃▆
income 0 1 6797.90 4245.92 611.00 4106.00 5930.50 8187.25 25879.00 ▇▇▁▁▁
women 0 1 28.98 31.72 0.00 3.59 13.60 52.20 97.51 ▇▂▂▂▂
prestige 0 1 46.83 17.20 14.80 35.23 43.60 59.28 87.20 ▃▇▅▅▂
census 0 1 5401.77 2644.99 1113.00 3120.50 5135.00 8312.50 9517.00 ▆▆▃▂▇

Letra A

  • Faça o gráfico de dispersão da variável resposta prestige pelas variáveis explicativas income e education.
ggplot2::ggplot(carData::Prestige, ggplot2::aes(x = income, y = prestige)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Renda média em Doláres") + 
  ggplot2::ylab("Score do prestígio") +
  ggplot2::labs(title = "Gráfico de dispersão entre o Score do Prestígio Pineo-Porter com relação a renda média", subtitle = "Pesquisa social realizada em meados dos anos 60.")

ggplot2::ggplot(carData::Prestige, ggplot2::aes(x = education, y = prestige)) +
  ggplot2::geom_point()+
  ggplot2::xlab("média em anos de estudo") + 
  ggplot2::ylab("Score do prestígio") +
  ggplot2::labs(title = "Gráfico de dispersão entre o Score do Prestígio Pineo-Porter com relação a média em anos de estudo", subtitle = "Pesquisa social realizada em meados dos anos 60.")

Letra B

  • Realize o ajuste de um modelo GAM com a variável resposta prestige tendo uma distribuição Normal. Faça o gráfico das funções de suavização.

Verificando a distribuição dos dados da variável prestige, se tem um bom ajuste considerando a distribuição normal.

ggplot2::ggplot(carData::Prestige, ggplot2::aes(prestige))+
  ggplot2::geom_histogram(fill='white', 
                          color = "black", 
                          breaks = hist(carData::Prestige$prestige, plot = F)$breaks) +
  ggplot2::xlab("Score do Prestígio de Pineo-Porter")+
  ggplot2::ylab("Frequência")+
  ggplot2::labs(title = "Histograma do Score de Prestígio de Pineo-Porter para a ocupação", subtitle = "Pesquisa social realizada em meados dos anos 60")

Vamos continuar o ajuste considerando uma distribuição normal.

Primeiro Ajuste

fit3 <- mgcv::gam(prestige~te(income)+te(education), family = gaussian(), data = carData::Prestige)
summary(fit3)

Family: gaussian 
Link function: identity 

Formula:
prestige ~ te(income) + te(education)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  46.8333     0.6858   68.29   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                edf Ref.df    F p-value    
te(income)    2.843  3.264 18.5  <2e-16 ***
te(education) 2.912  3.412 44.9  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.838   Deviance explained = 84.7%
GCV = 51.379  Scale est. = 47.976    n = 102
par(mfrow=c(1,2))
plot(fit3)

Segundo ajuste

fit4 <- mgcv::gam(prestige~t2(income)+t2(education), family = gaussian(), data = carData::Prestige)
summary(fit4)

Family: gaussian 
Link function: identity 

Formula:
prestige ~ t2(income) + t2(education)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  46.8333     0.6858   68.29   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                edf Ref.df      F  p-value    
t2(income)    2.843  3.264  7.189 9.05e-05 ***
t2(education) 2.912  3.412 44.249  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.838   Deviance explained = 84.7%
GCV = 51.379  Scale est. = 47.976    n = 102
par(mfrow=c(1,2))
plot(fit4)

Terceiro Ajuste

fit5 <- mgcv::gam(prestige~s(income)+s(education), family = gaussian(), data = carData::Prestige)
summary(fit5)

Family: gaussian 
Link function: identity 

Formula:
prestige ~ s(income) + s(education)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  46.8333     0.6889   67.98   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
               edf Ref.df     F p-value    
s(income)    3.118  3.877 14.61  <2e-16 ***
s(education) 3.177  3.952 38.78  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.836   Deviance explained = 84.7%
GCV = 52.143  Scale est. = 48.414    n = 102
par(mfrow=c(1,2))
plot(fit5)

Resultado: Conforme os ajustes acima, o segundo e primeiro ajuste, no qual utiliza o parâmetro de suavização t2 e t, se mostraram mais eficiente, tendo um coeficiente R^2 ajustado de 83.8% e um Desvio explicado de cerca de 84.7%, sendo assim utilizaremos o segundo modelo ajustado.

Letra C

  • Faça uma análise de diagnósticos do modelo escolhido. O que você pode concluir do modelo?
par(mfrow = c(2,2))
mgcv::gam.check(fit4)


Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 4 iterations.
The RMS GCV score gradient at convergence was 0.000514641 .
The Hessian was positive definite.
Model rank =  9 / 9 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                k'  edf k-index p-value
t2(income)    4.00 2.84    0.99    0.40
t2(education) 4.00 2.91    1.02    0.55

Resposta: De acordo com a análise gráfica dos resíduos podemos verificar que o modelo ficou bem ajustado e que os desvios são aproximadamente normais.

Questão 3

  • Considere o banco de dados fabric do pacote gamlss do R. Em que y é o número de falhas em um rolo de tecido e leng é o comprimento do tecido. A variável x, que é o log de leng não usaremos na questão.

Verificando a categoria das variáveis no banco

glimpse(gamlss.data::fabric)
Rows: 32
Columns: 3
$ leng <dbl> 551, 651, 832, 375, 715, 868, 271, 630, 491, 372, 645, 441, 895, …
$ y    <dbl> 6, 4, 17, 9, 14, 8, 5, 7, 7, 7, 6, 8, 28, 4, 10, 4, 8, 9, 23, 9, …
$ x    <dbl> 6.311735, 6.478510, 6.723832, 5.926926, 6.572283, 6.766192, 5.602…

Estatística Descritiva

Visualizando as primeiras linhas do banco

gamlss.data::fabric %>%
  slice_head(n=10) %>%
  gt::gt()
leng y x
551 6 6.311735
651 4 6.478510
832 17 6.723832
375 9 5.926926
715 14 6.572283
868 8 6.766192
271 5 5.602119
630 7 6.445720
491 7 6.196444
372 7 5.918894
 
skimr::skim(gamlss.data::fabric)
Data summary
Name gamlss.data::fabric
Number of rows 32
Number of columns 3
_______________________
Column type frequency:
numeric 3
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
leng 0 1 587.66 211.39 122.0 453.75 590.50 735.75 952.00 ▂▃▆▇▅
y 0 1 8.88 5.81 1.0 5.75 8.00 9.25 28.00 ▆▇▂▁▁
x 0 1 6.29 0.46 4.8 6.12 6.38 6.60 6.86 ▁▁▂▅▇

Letra A

  • Faça o gráfico de dispersão da variável resposta y pela variável explicativa (x).
ggplot2::ggplot(gamlss.data::fabric, ggplot2::aes(x = leng, y = y)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Comprimento do tecido") + 
  ggplot2::ylab("Número de falhas") +
  ggplot2::labs(title = "Gráfico de dispersão entre o número de falhas em um rolo e o comprimento do tecido")

Letra B

  • Realize o ajuste de um modelo GAMLSS com a variável resposta y tendo uma distribuição Poisson.

Verificando a distribuição dos dados da variável y, se tem um bom ajuste considerando a distribuição Poisson.

ggplot2::ggplot(gamlss.data::fabric, ggplot2::aes(y))+
  ggplot2::geom_histogram(fill='white', 
                          color = "black", 
                          breaks = hist(gamlss.data::fabric$y, plot = F)$breaks) +
  ggplot2::xlab("Número de falhas")+
  ggplot2::ylab("Frequência")+
  ggplot2::labs(title = "Histograma do número de falhas em um rolo de tecido")

Vamos ajustar o modelo de acordo com a distribuição Poisson para a variável resposta.

fit6 <- gamlss::gamlss(y~leng, sigma.formula = ~leng, data = gamlss.data::fabric, family = PO)
GAMLSS-RS iteration 1: Global Deviance = 185.0559 
GAMLSS-RS iteration 2: Global Deviance = 185.0559 
summary(fit6)
Warning in summary.gamlss(fit6): summary: vcov has failed, option qr is used instead
******************************************************************
Family:  c("PO", "Poisson") 

Call:  gamlss::gamlss(formula = y ~ leng, sigma.formula = ~leng, family = PO,  
    data = gamlss.data::fabric) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  log
Mu Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 0.9717506  0.2124710   4.574 4.79e-06 ***
leng        0.0019297  0.0003063   6.300 2.97e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  32 
Degrees of Freedom for the fit:  2
      Residual Deg. of Freedom:  30 
                      at cycle:  2 
 
Global Deviance:     185.0559 
            AIC:     189.0559 
            SBC:     191.9873 
******************************************************************

Letra C

  • Faça uma análise de diagnósticos do modelo escolhido. O que você pode concluir do modelo?
plot(fit6)

******************************************************************
     Summary of the Randomised Quantile Residuals
                           mean   =  -0.04398648 
                       variance   =  1.991556 
               coef. of skewness  =  0.4198306 
               coef. of kurtosis  =  2.800601 
Filliben correlation coefficient  =  0.9893308 
******************************************************************

Resposta: Considerando que o problema tem poucas observações e que isto não traz muita precisão ao ajuste do modelo, todavia, concluimos que o modelo foi bem ajustado.