Regressão 2 - Aula 8 e 9

Author

Paulo Manoel da Silva Junior

Aula 8 - Uso do GAMLS

Os modelos aditivos generalizados para locação, escala e forma (GAMLSS – Generalized Additive Models for Location, Scale and Shape) foram introduzidos por Rigby e Stasinopoulos (2001, 2005) para compensar as limitações dos MLGs e GAMs.

Banco de Dados

  • Em nosso exemplo, utilizaremos o banco de dados usair Poluição do ar dos EUA retirados de Hand et al (1994). Contém as seguintes variáveis:

y: concentração de dióxido de enxofre no ar em mgs por metro cubico em 41 cidades;

x1: temperatura média anual em graus F;

x2: número de fabricantes que empregam mais que 20 trabalhadores;

x3: tamanho da população em milhares;

x4: velocidade média anual do vento em milhas por hora;

x5: precipitação média anual em polegadas;

x6: número médio de dias de chuva por ano.

Estatístitca Descritiva

Visualizando algumas linhas do banco de dados

gamlss.data::usair %>%
  slice_head(n=15) %>%
  gt::gt()
y x1 x2 x3 x4 x5 x6
10 70.3 213 582 6.0 7.05 36
13 61.0 91 132 8.2 48.52 100
12 56.7 453 716 8.7 20.66 67
17 51.9 454 515 9.0 12.95 86
56 49.1 412 158 9.0 43.37 127
36 54.0 80 80 9.0 40.25 114
29 57.3 434 757 9.3 38.89 111
14 68.4 136 529 8.8 54.47 116
10 75.5 207 335 9.0 59.80 128
24 61.5 368 497 9.1 48.34 115
110 50.6 3344 3369 10.4 34.44 122
28 52.3 361 746 9.7 38.74 121
17 49.0 104 201 11.2 30.85 103
8 56.6 125 277 12.7 30.58 82
30 55.6 291 593 8.3 43.11 123

Verificando o tipo das variáveis

glimpse(gamlss.data::usair)
Rows: 41
Columns: 7
$ y  <int> 10, 13, 12, 17, 56, 36, 29, 14, 10, 24, 110, 28, 17, 8, 30, 9, 47, …
$ x1 <dbl> 70.3, 61.0, 56.7, 51.9, 49.1, 54.0, 57.3, 68.4, 75.5, 61.5, 50.6, 5…
$ x2 <int> 213, 91, 453, 454, 412, 80, 434, 136, 207, 368, 3344, 361, 104, 125…
$ x3 <int> 582, 132, 716, 515, 158, 80, 757, 529, 335, 497, 3369, 746, 201, 27…
$ x4 <dbl> 6.0, 8.2, 8.7, 9.0, 9.0, 9.0, 9.3, 8.8, 9.0, 9.1, 10.4, 9.7, 11.2, …
$ x5 <dbl> 7.05, 48.52, 20.66, 12.95, 43.37, 40.25, 38.89, 54.47, 59.80, 48.34…
$ x6 <int> 36, 100, 67, 86, 127, 114, 111, 116, 128, 115, 122, 121, 103, 82, 1…
visdat::vis_dat(gamlss.data::usair)

Observando as correlações

visdat::vis_cor(gamlss.data::usair)

Medidas Descritivas do banco de dados

skimr::skim(gamlss.data::usair)
Data summary
Name gamlss.data::usair
Number of rows 41
Number of columns 7
_______________________
Column type frequency:
numeric 7
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
y 0 1 30.05 23.47 8.00 13.00 26.00 35.00 110.0 ▇▃▂▁▁
x1 0 1 55.76 7.23 43.50 50.60 54.60 59.30 75.5 ▃▇▅▂▁
x2 0 1 463.10 563.47 35.00 181.00 347.00 462.00 3344.0 ▇▁▁▁▁
x3 0 1 608.61 579.11 71.00 299.00 515.00 717.00 3369.0 ▇▂▁▁▁
x4 0 1 9.44 1.43 6.00 8.70 9.30 10.60 12.7 ▁▂▇▅▁
x5 0 1 36.77 11.77 7.05 30.96 38.74 43.11 59.8 ▂▁▆▇▁
x6 0 1 113.90 26.51 36.00 103.00 115.00 128.00 166.0 ▁▂▅▇▂

Verificando qual a distribução dos dados através do histograma

A variável dependente é a y, que se trata de sulpher dioxide concentration in air mgs. per cubic metre in 41 cities in the USA

ggplot2::ggplot(gamlss.data::usair, ggplot2::aes(x = y))+
  ggplot2::geom_histogram(fill='white', 
                          color = "black", 
                          breaks = hist(gamlss.data::usair$y, plot = F)$breaks)  + ggplot2::xlab("sulpher dioxide concentration in air mgs") + ggplot2::ylab("Frequência") + ggplot2::labs(title = "Sulpher dioxide concentration in air mgs", subtitle = "per cubic metre in 41 cities in the USA")

Como é observado vamos utilizar a distribuição Gamma, a qual se adapta melhor aos dados.

Modelagem

Seleção de Variáveis

Seleção de Variáveis para Mu

fit1_mu <- gamlss::gamlss(y~1, sigma.formula = ~1, data = gamlss.data::usair, family = GA)
GAMLSS-RS iteration 1: Global Deviance = 349.7146 
GAMLSS-RS iteration 2: Global Deviance = 349.7146 
step_mu <- gamlss::stepGAIC(fit1_mu, scope = list(lower = ~1, upper =~ x1+x2+x3+x4+x5+x6))
Distribution parameter:  mu 
Start:  AIC= 353.71 
 y ~ 1 

       Df    AIC
+ x1    1 338.04
+ x2    1 343.05
+ x6    1 343.97
+ x3    1 349.20
<none>    353.71
+ x4    1 355.02
+ x5    1 355.43

Step:  AIC= 338.04 
 y ~ x1 

       Df    AIC
+ x2    1 328.78
+ x3    1 332.94
+ x5    1 335.11
+ x6    1 335.99
<none>    338.04
+ x4    1 339.65
- x1    1 353.71

Step:  AIC= 328.78 
 y ~ x1 + x2 

       Df    AIC
+ x5    1 324.44
+ x6    1 325.54
+ x3    1 325.67
+ x4    1 327.65
<none>    328.78
- x2    1 338.04
- x1    1 343.05

Step:  AIC= 324.44 
 y ~ x1 + x2 + x5 

       Df    AIC
+ x4    1 319.39
+ x3    1 322.48
<none>    324.44
+ x6    1 326.16
- x5    1 328.78
- x2    1 335.11
- x1    1 344.28

Step:  AIC= 319.39 
 y ~ x1 + x2 + x5 + x4 

       Df    AIC
+ x3    1 317.16
<none>    319.39
+ x6    1 321.39
- x4    1 324.44
- x5    1 327.65
- x2    1 335.43
- x1    1 346.16

Step:  AIC= 317.16 
 y ~ x1 + x2 + x5 + x4 + x3 

       Df    AIC
<none>    317.16
+ x6    1 319.16
- x3    1 319.39
- x4    1 322.48
- x5    1 324.14
- x2    1 324.92
- x1    1 336.11

Podemos observar que as variáveis que foram selecionadas para \mu, foram x1, x2, x3, x4 e x5.

Seleção de Variáveis para Sigma

fit2_sigma <- gamlss::gamlss(y~x1+x2+x3+x4+x5, data = gamlss.data::usair, family = GA) 
GAMLSS-RS iteration 1: Global Deviance = 303.162 
GAMLSS-RS iteration 2: Global Deviance = 303.1619 
step_sigma <- gamlss::stepGAIC(fit2_sigma, what = "sigma", scope = list(lower=~1, upper = ~ x1+x2+x3+x4+x5))
Distribution parameter:  sigma 
Start:  AIC= 317.16 
 ~1 

       Df    AIC
+ x5    1 315.58
+ x4    1 315.96
<none>    317.16
+ x3    1 317.39
+ x1    1 317.79
+ x2    1 318.61

Step:  AIC= 315.58 
 ~x5 
Warning in RS(): Algorithm RS has not yet converged
       Df    AIC
+ x1    1 314.08
<none>    315.58
+ x3    1 315.63
+ x2    1 317.07
- x5    1 317.16
+ x4    1 317.77

Step:  AIC= 314.08 
 ~x5 + x1 

       Df    AIC
<none>    314.08
+ x3    1 314.48
+ x2    1 315.38
+ x4    1 315.51
- x1    1 315.58
- x5    1 317.79

Podemos observar que as variáveis que foram selecionadas para \sigma foram as variáveis x1 e x5.

Ajuste Final

Procedendo com o ajuste final, já que a seleção de variáveis tanto para \mu, como para \sigma foram bem ajustadas.

fit_final <- gamlss::gamlss(y~x1+x2+x3+x4+x5, sigma.formula = ~ x1+x5, data = gamlss.data::usair, family = GA)
GAMLSS-RS iteration 1: Global Deviance = 300.755 
GAMLSS-RS iteration 2: Global Deviance = 298.3819 
GAMLSS-RS iteration 3: Global Deviance = 297.0355 
GAMLSS-RS iteration 4: Global Deviance = 296.3957 
GAMLSS-RS iteration 5: Global Deviance = 296.1655 
GAMLSS-RS iteration 6: Global Deviance = 296.1034 
GAMLSS-RS iteration 7: Global Deviance = 296.0886 
GAMLSS-RS iteration 8: Global Deviance = 296.0852 
GAMLSS-RS iteration 9: Global Deviance = 296.0845 
summary(fit_final)
Warning in summary.gamlss(fit_final): summary: vcov has failed, option qr is used instead
******************************************************************
Family:  c("GA", "Gamma") 

Call:  
gamlss::gamlss(formula = y ~ x1 + x2 + x3 + x4 + x5, sigma.formula = ~x1 +  
    x5, family = GA, data = gamlss.data::usair) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  log
Mu Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.5677204  0.7072040  10.701 1.40e-12 ***
x1          -0.0557872  0.0088892  -6.276 3.35e-07 ***
x2           0.0013030  0.0004403   2.959   0.0055 ** 
x3          -0.0007487  0.0004183  -1.790   0.0821 .  
x4          -0.2244998  0.0378302  -5.934 9.42e-07 ***
x5           0.0194673  0.0035151   5.538 3.14e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.28178    0.84700   0.333  0.74120    
x1          -0.04962    0.01642  -3.022  0.00447 ** 
x5           0.04064    0.01018   3.991  0.00029 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  41 
Degrees of Freedom for the fit:  9
      Residual Deg. of Freedom:  32 
                      at cycle:  9 
 
Global Deviance:     296.0845 
            AIC:     314.0845 
            SBC:     329.5066 
******************************************************************

Análise gráfica

plot(fit_final)

******************************************************************
          Summary of the Quantile Residuals
                           mean   =  0.04304783 
                       variance   =  1.051775 
               coef. of skewness  =  0.2085559 
               coef. of kurtosis  =  2.420059 
Filliben correlation coefficient  =  0.9951659 
******************************************************************

Conclusão

Podemos observar através do resultado do ajuste que foi fornecido que o modelo foi bem ajustado, e a análise gráfica comprova, pois os resíduos aparentam estar bem distribuidos, seguindo uma distribuição normal, o que é evidenciado pela densidade que foi plotada.

Aula 9 - Regressão Não Linear

Os modelos de regressão linear fornecem uma estrutura rica e flexível que supre as necessidades na solução de vários problemas práticos. Entretanto, os modelos de regressão linear não são apropriados em todas as situações.

Por exemplo, em problemas nas áreas da engenharia e das ciências, a relação entre a variável resposta e as variáveis explicativas está relacionada com uma função não-linear conhecida. Além disso, os modelos não-lineares são deduzidos a partir de suposições teóricas, sendo os parâmetros resultantes interpretáveis de modo que aproximá-los por um modelo linear prejudicaria bastante a obtenção de estimativas mais realistas dos parâmetros de interesse.

Banco de Dados

O banco de dados fornece o peso, em quilogramas, de um paciente obeso em 52 pontos de tempo durante um período de 8 meses de um programa de reabilitação de peso.

Observando algumas linhas do banco de dados

MASS::wtloss %>%
  slice_head(n=15) %>%
  gt::gt()
Days Weight
0 184.35
4 182.51
7 180.45
7 179.91
11 177.91
18 175.81
24 173.11
30 170.06
32 169.31
43 165.10
46 163.11
60 158.30
64 155.80
70 154.31
71 153.86

Estatístitca Descritiva

skimr::skim(MASS::wtloss)
Data summary
Name MASS::wtloss
Number of rows 52
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
Days 0 1 119.15 71.52 0.0 68.5 127.5 165.00 246.00 ▇▇▇▇▆
Weight 0 1 142.25 21.10 112.6 127.5 135.8 154.68 184.35 ▇▇▆▂▅

Visualização Gráfica

ggplot2::ggplot(MASS::wtloss, ggplot2::aes(x = Days, y = Weight)) +
  ggplot2::geom_point()+
  ggplot2::geom_smooth(method = "gam") +
  ggplot2::xlab("Dia") + ggplot2::ylab("Peso") + 
  ggplot2::labs(title = "Peso de um paciente obeso durante 8 meses", subtitle = "Em um período de reabilitação de peso")
`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

Ajustando o modelo

  • Deve ser lembrando que precisamos saber da função e do chute inicial para que a modelagem possa acontecer
wtloss.fm <- nls(Weight ~ b0 + b1*2^(-Days/th), data = MASS::wtloss,start = list(b0=90, b1=95, th=120))
summary(wtloss.fm)

Formula: Weight ~ b0 + b1 * 2^(-Days/th)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b0   81.374      2.269   35.86   <2e-16 ***
b1  102.684      2.083   49.30   <2e-16 ***
th  141.911      5.295   26.80   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8949 on 49 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 4.324e-06

Análise gráfica do ajuste

plot(MASS::wtloss)
with(wtloss, lines(Days, fitted(wtloss.fm)))
abline(wtloss.fm,lty = 2)
Warning in abline(wtloss.fm, lty = 2): usando somente os dois primeiros de 3
coeficientes de regressão