MLGExemplo03

O Cardiovascular Health Study (CHS) é um estudo longitudinal de base populacional de doença coronariana e derrame cerebral em adultos com 65 anos ou mais. Os participantes do estudo foram recrutados em 1989-1990 em quatro comunidades: Forsyth County, NC; Condado de Sacramento, CA; Condado de Washington, MD; e Pittsburgh, PA. Os participantes elegíveis foram retirados das listas de elegibilidade do Medicare em cada área. Pensa-se que essas listas estão 98-99% completas para cidadãos dos EUA com 65 anos ou mais. Os participantes elegíveis incluíram todas as pessoas que viviam na casa de cada indivíduo da amostra das listas que tinham: 1) 65 anos ou mais; 2) não institucionalizado; 3) expectativa de permanência na área por 3 anos; e 4) capaz de dar consentimento informado e não exigiu um respondente por procuração. Os participantes eram elegíveis para CHS quando tinham ou não história de doença cardiovascular. No entanto, as pessoas foram excluídas se estivessem em cadeiras de rodas ou se estivessem recebendo tratamento hospitalar ou de câncer.

Os dados para este estudo consistem no subconjunto de participantes que foram recrutados na primeira onda de recrutamento e eram ‘saudáveis’, ou seja, não tinham histórico de doenças cardíacas ou circulatórias, nenhuma restrição de atividades diárias por doença e nenhum medicamento que indicasse doença cardíaca. Um grande número de variáveis, foi determinado para cada participante do estudo na linha de base, ou seja, no momento do recrutamento. O exame inicial consistiu em uma entrevista domiciliar e exame clínico. Durante a entrevista domiciliar, foram coletadas informações sobre histórico médico anterior, uso médico e atividade física. Também foram obtidas informações sobre a presença de alterações do funcionamento físico. O exame clínico incluiu coleta de sangue em jejum e medições da pressão arterial em posição sentada. A tabela abaixo contém uma descrição das variáveis, no conjunto de dados que você usará. Esses dados podem ser encontrados em https://www.ics.uci.edu/~dgillen/STAT211/Data/chsData.txt.

Variável		Descrição
clinic		1 = Sacremento, 2 = Forsyth, 3 = Washington, 4 = Pittsburgh
initdate		Data de recrutamento (em dias, medida a partir de uma data de início arbitrária)
season		Estação inicial, 1 = verão, 2 = outono, 3 = inverno, 4 = primavera
gender		0 = Feminino, 1 = Masculino
age		Idade na linha de base (em anos)
weight		Peso na linha de base (em libras)
weight50		Peso recuperado aos 50 anos (em lbs)
grade		Anos de educação
arth		1 = artrite, 0 = nenhuma (na linha de base)
sbp		Pressão arterial sistólica basal (em mmHg)
pkyrs		Maço de anos de história de tabagismo (número de anos como fumante × número de maços / dia × 365
diab		Diabetes no início do estudo, 1 = nenhum, 2 = limítrofe, 3 = diabetes
income		Renda familiar (em 1k dólares), 1 = ≤ 5, 2 = 5-8, 3 = 8-12, 4 = 12-16, 5 = 16-24, 6 = 24-35, 7 = 35-50, 8 = ≥ 50
exint0		Medida basal da intensidade do exercício, 0 = nenhum exercício, 1 = baixa intensidade, 2 = intensidade moderada, 3 = intensidade alta
block0		Medida de linha de base de blocos percorridos nas últimas 2 semanas (em cerca de 12 blocos/milha)
kcal0		Medida de linha de base de quilocalorias estimadas gastas na atividade de exercício nas últimas 2 semanas

Dados

O objetivo desta análise é investigar a associação entre a resposta, pressão arterial sistólica maior que 140 mmHg e o preditor de peso de interesse. Em particular, há a hipótese de que o aumento do peso está associado ao aumento da pressão arterial.

dados = read.table("https://www.ics.uci.edu/~dgillen/STAT211/Data/chsData.txt", header = T)
head(dados)

##   clinic initdate season gender age weight weight50 grade arth sbp pkyrs diab
## 1      1   148665      2      1  76  213.5      212     9    1 123    27    1
## 2      1   148772      3      0  67  140.0      130    11    1 159     0    1
## 3      1   148706      2      0  79  140.0      155    20    0 135     0    1
## 4      1   148701      2      1  86  135.0      138    21    0 142     0    1
## 5      1   148859      4      1  70  193.0      215    20    0 128    33    2
## 6      1   148604      1      0  68  153.0      135    20    0 118     0    1
##   income exint0 block0  kcal0
## 1      5      2     60 1005.0
## 2      4      0      0    0.0
## 3      6      1      6  240.0
## 4      6      1     30  420.0
## 5      6      2     15 2977.5
## 6      7      1     24  210.0

library(lattice)
xyplot(sbp ~ age | season, groups = gender, data = dados,
       auto.key=list(space="top", columns=2, title="Sexo", cex.title=1))

xyplot(sbp ~ age | season, groups = exint0, data = dados,
       auto.key=list(space="top", columns=4, title="Intensidade do exercício", cex.title=1))

xyplot(sbp>140 ~ age | season, groups = gender, data = dados,
       auto.key=list(space="top", columns=2, title="Sexo", cex.title=1))

xyplot(sbp>140 ~ pkyrs | season, groups = season, data = dados,
       auto.key=list(space="top", columns=4, title="Estação", cex.title=1))

xyplot(sbp>140 ~ weight | season, groups = season, data = dados,
       auto.key=list(space="top", columns=4, title="Estação", cex.title=1))

xyplot(sbp>140 ~ weight50 | season, groups = season, data = dados,
       auto.key=list(space="top", columns=4, title="Estação", cex.title=1))

Não apresenta-se nehuma variável claramente como explicativa isolada do sbp, a pressão arterial sistólica basal. Vejamos o que acontece quando mostramos a pressão arterial associada ao aumento de peso.

xyplot(sbp>140 ~ weight | season, groups = exint0, data = dados,
       auto.key=list(space="top", columns=4, title="Intensidade do exercício", cex.title=1))

xyplot(sbp ~ weight50 | season, groups = exint0, data = dados,
       auto.key=list(space="top", columns=4, title="Intensidade do exercício", cex.title=1))

Modelos

Primeiro pensemos na resposta sbp como contínua, iso significa que podemos atribuir à resposta distribuições Normal, Gama, etc. Vejamos estas duas.

library(gamlss)

Temos uma situação nesta base de dados: muitas das variáveis explicativas estão codificadas numéricamente, isso significa que se não indicamos que são fatores, serão consideradas como numéricas na estimação.

Faremos esse trabalho agora, desconsiderando algumas que imaginamos sejam não informativas: clinic e season,

Modelo Normal

modelo0 = gamlss(sbp ~ initdate + factor(gender) + age + weight + weight50 + grade + factor(arth) +
                   pkyrs + factor(diab) + factor(income) + factor(exint0) + block0 + kcal0, 
                 data = na.omit(dados), family = NO(),
                 control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(modelo0)

## Warning in summary.gamlss(modelo0): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = sbp ~ initdate + factor(gender) + age + weight +  
##     weight50 + grade + factor(arth) + pkyrs + factor(diab) +  
##     factor(income) + factor(exint0) + block0 + kcal0, family = NO(),  
##     data = na.omit(dados), control = gamlss.control(n.cyc = 2000,  
##         trace = FALSE)) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.192e+03  6.769e+02  -3.239  0.00122 ** 
## initdate         1.521e-02  4.558e-03   3.338  0.00086 ***
## factor(gender)1 -5.647e-02  1.254e+00  -0.045  0.96410    
## age              8.288e-01  9.650e-02   8.588  < 2e-16 ***
## weight           1.146e-01  2.763e-02   4.147 3.51e-05 ***
## weight50        -8.796e-02  3.232e-02  -2.722  0.00655 ** 
## grade            2.748e-02  1.140e-01   0.241  0.80951    
## factor(arth)1    1.460e+00  9.115e-01   1.602  0.10933    
## pkyrs           -1.792e-02  1.827e-02  -0.981  0.32689    
## factor(diab)2    3.590e+00  1.373e+00   2.615  0.00899 ** 
## factor(diab)3    7.255e+00  1.504e+00   4.824 1.51e-06 ***
## factor(income)2 -3.101e+00  3.376e+00  -0.919  0.35846    
## factor(income)3 -2.356e+00  3.196e+00  -0.737  0.46115    
## factor(income)4 -1.944e-01  3.098e+00  -0.063  0.94997    
## factor(income)5 -1.354e+00  3.052e+00  -0.444  0.65742    
## factor(income)6 -2.570e+00  3.108e+00  -0.827  0.40839    
## factor(income)7 -1.403e+00  3.240e+00  -0.433  0.66516    
## factor(income)8 -2.886e+00  3.223e+00  -0.895  0.37065    
## factor(exint0)1  1.187e+00  2.119e+00   0.560  0.57538    
## factor(exint0)2 -1.394e+00  2.188e+00  -0.637  0.52408    
## factor(exint0)3 -1.218e+00  2.492e+00  -0.489  0.62506    
## block0          -2.259e-03  8.464e-03  -0.267  0.78962    
## kcal0            4.743e-04  3.027e-04   1.567  0.11729    
## ---
## 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)  3.00902    0.01557   193.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  2063 
## Degrees of Freedom for the fit:  24
##       Residual Deg. of Freedom:  2039 
##                       at cycle:  2 
##  
## Global Deviance:     18269.77 
##             AIC:     18317.77 
##             SBC:     18452.94 
## ******************************************************************

plot(modelo0)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -4.484391e-16 
##                        variance   =  1.000485 
##                coef. of skewness  =  0.5326837 
##                coef. of kurtosis  =  3.351368 
## Filliben correlation coefficient  =  0.9916229 
## ******************************************************************

Modelo Gama

modelo1 = update(modelo0, family = GA())
summary(modelo1)

## Warning in summary.gamlss(modelo1): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = sbp ~ initdate + factor(gender) + age + weight +  
##     weight50 + grade + factor(arth) + pkyrs + factor(diab) +  
##     factor(income) + factor(exint0) + block0 + kcal0, family = GA(),  
##     data = na.omit(dados), control = gamlss.control(n.cyc = 2000,  
##         trace = FALSE)) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -1.227e+01  4.950e+00  -2.479 0.013273 *  
## initdate         1.122e-04  3.333e-05   3.367 0.000774 ***
## factor(gender)1 -4.688e-04  9.174e-03  -0.051 0.959247    
## age              6.128e-03  7.057e-04   8.683  < 2e-16 ***
## weight           8.565e-04  2.021e-04   4.239 2.35e-05 ***
## weight50        -6.490e-04  2.363e-04  -2.746 0.006078 ** 
## grade            2.070e-04  8.336e-04   0.248 0.803915    
## factor(arth)1    1.103e-02  6.666e-03   1.655 0.098170 .  
## pkyrs           -1.313e-04  1.336e-04  -0.983 0.325927    
## factor(diab)2    2.702e-02  1.004e-02   2.691 0.007179 ** 
## factor(diab)3    5.252e-02  1.100e-02   4.775 1.92e-06 ***
## factor(income)2 -2.291e-02  2.469e-02  -0.928 0.353420    
## factor(income)3 -1.726e-02  2.337e-02  -0.739 0.460193    
## factor(income)4 -1.557e-03  2.266e-02  -0.069 0.945217    
## factor(income)5 -1.048e-02  2.232e-02  -0.469 0.638773    
## factor(income)6 -1.867e-02  2.273e-02  -0.822 0.411429    
## factor(income)7 -1.053e-02  2.370e-02  -0.445 0.656718    
## factor(income)8 -2.089e-02  2.357e-02  -0.886 0.375550    
## factor(exint0)1  9.412e-03  1.549e-02   0.607 0.543605    
## factor(exint0)2 -9.964e-03  1.600e-02  -0.623 0.533440    
## factor(exint0)3 -9.851e-03  1.822e-02  -0.541 0.588830    
## block0          -2.139e-05  6.189e-05  -0.346 0.729640    
## kcal0            3.636e-06  2.214e-06   1.642 0.100663    
## ---
## 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) -1.90909    0.01551  -123.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  2063 
## Degrees of Freedom for the fit:  24
##       Residual Deg. of Freedom:  2039 
##                       at cycle:  2 
##  
## Global Deviance:     18175.98 
##             AIC:     18223.98 
##             SBC:     18359.15 
## ******************************************************************

plot(modelo1)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.0001508886 
##                        variance   =  1.000469 
##                coef. of skewness  =  0.2451554 
##                coef. of kurtosis  =  2.947309 
## Filliben correlation coefficient  =  0.9977928 
## ******************************************************************

Perece-nos mais adequado, visualmente, o modelo gama mesmo com o erro no resumo do modelo ajustado. O padrão na comando summary.gamlss é type = “vcov” usa utiliza o método vcov() para gamlss para obter a matriz de variância-covariância dos coeficientes da regressão estimados. A alternativa type = “qr” é o método original usado em gamlss para estimar os erros padrão, mas não é confiável, pois não leva em consideração a inter-correlação entre os parâmetros de distribuição.

A matriz de variância e covariância é calculada usando o inverso das segundas derivadas numéricas da matriz de informação observada. Este é um método mais confiável, pois leva em consideração a inter-correlação entre todos os parâmetros. O type = “qr” assume que os parâmetros são fixos nos valores estimados. Observe que ambos os métodos não são apropriados e devem ser usados com cuidado se os termos de suavização forem usados no ajuste.

Modelo Binomial

Mudamos agora nosso pensamento. Como sugerido na descrição dos dados, queremos investigar a associação entre a resposta, pressão arterial sistólica maior que 140 mmHg e o preditor de peso de interesse. Em particular, há a hipótese de que o aumento do peso está associado ao aumento da pressão arterial.

Para isto faremos a codificação da reposta como pressão arterial sistólica maior do que 140 mmHg ou não, segundo todos os outros preditores considerados até agora.

modelo2 = gamlss(sbp>140 ~ initdate + factor(gender) + age + weight + weight50 + grade + factor(arth) +
                   pkyrs + factor(diab) + factor(income) + factor(exint0) + block0 + kcal0, 
                 data = na.omit(dados), family = BI(),
                 control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(modelo2)

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = sbp > 140 ~ initdate + factor(gender) +  
##     age + weight + weight50 + grade + factor(arth) +  
##     pkyrs + factor(diab) + factor(income) + factor(exint0) +  
##     block0 + kcal0, family = BI(), data = na.omit(dados),  
##     control = gamlss.control(n.cyc = 2000, trace = FALSE)) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)     -1.876e+02  9.553e-01 -196.331  < 2e-16 ***
## initdate         1.226e-03  7.465e-06  164.179  < 2e-16 ***
## factor(gender)1 -2.720e-03  1.307e-01   -0.021  0.98340    
## age              6.509e-02  1.007e-02    6.461 1.30e-10 ***
## weight           8.503e-03  2.903e-03    2.929  0.00344 ** 
## weight50        -8.084e-03  3.411e-03   -2.370  0.01787 *  
## grade            1.608e-04  1.198e-02    0.013  0.98929    
## factor(arth)1    1.542e-01  9.600e-02    1.607  0.10828    
## pkyrs           -2.899e-04  1.929e-03   -0.150  0.88056    
## factor(diab)2    2.872e-01  1.419e-01    2.024  0.04315 *  
## factor(diab)3    6.062e-01  1.528e-01    3.968 7.51e-05 ***
## factor(income)2 -2.377e-01  3.487e-01   -0.682  0.49557    
## factor(income)3 -1.410e-01  3.299e-01   -0.427  0.66924    
## factor(income)4  5.481e-02  3.188e-01    0.172  0.86352    
## factor(income)5  5.244e-02  3.142e-01    0.167  0.86747    
## factor(income)6 -2.187e-01  3.215e-01   -0.680  0.49638    
## factor(income)7 -1.192e-01  3.354e-01   -0.356  0.72221    
## factor(income)8 -2.370e-01  3.342e-01   -0.709  0.47836    
## factor(exint0)1 -1.584e-01  2.173e-01   -0.729  0.46621    
## factor(exint0)2 -3.089e-01  2.250e-01   -1.373  0.16997    
## factor(exint0)3 -2.528e-01  2.580e-01   -0.980  0.32743    
## block0          -6.893e-04  8.970e-04   -0.768  0.44231    
## kcal0            3.997e-05  4.680e-05    0.854  0.39321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  2063 
## Degrees of Freedom for the fit:  23
##       Residual Deg. of Freedom:  2040 
##                       at cycle:  2 
##  
## Global Deviance:     2587.528 
##             AIC:     2633.528 
##             SBC:     2763.062 
## ******************************************************************

plot(modelo2)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -9.13295e-05 
##                        variance   =  0.996911 
##                coef. of skewness  =  0.01801465 
##                coef. of kurtosis  =  2.967419 
## Filliben correlation coefficient  =  0.9997591 
## ******************************************************************

Modifiquemos a função de ligação, isso traz melhorias?

modelo3 = update(modelo2, family = BI(mu.link = "probit"))
summary(modelo3)

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = sbp > 140 ~ initdate + factor(gender) +  
##     age + weight + weight50 + grade + factor(arth) +  
##     pkyrs + factor(diab) + factor(income) + factor(exint0) +  
##     block0 + kcal0, family = BI(mu.link = "probit"),  
##     data = na.omit(dados), control = gamlss.control(n.cyc = 2000,  
##         trace = FALSE)) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  probit
## Mu Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)     -1.135e+02  5.817e-01 -195.154  < 2e-16 ***
## initdate         7.414e-04  7.465e-06   99.312  < 2e-16 ***
## factor(gender)1 -3.440e-03  8.061e-02   -0.043  0.96596    
## age              4.007e-02  6.120e-03    6.548 7.36e-11 ***
## weight           5.252e-03  1.770e-03    2.966  0.00305 ** 
## weight50        -4.957e-03  2.076e-03   -2.387  0.01706 *  
## grade            9.257e-05  7.305e-03    0.013  0.98989    
## factor(arth)1    9.576e-02  5.857e-02    1.635  0.10220    
## pkyrs           -1.828e-04  1.183e-03   -0.154  0.87726    
## factor(diab)2    1.774e-01  8.730e-02    2.032  0.04229 *  
## factor(diab)3    3.738e-01  9.428e-02    3.965 7.60e-05 ***
## factor(income)2 -1.508e-01  2.125e-01   -0.710  0.47793    
## factor(income)3 -8.779e-02  2.009e-01   -0.437  0.66218    
## factor(income)4  3.068e-02  1.942e-01    0.158  0.87449    
## factor(income)5  2.979e-02  1.914e-01    0.156  0.87632    
## factor(income)6 -1.355e-01  1.954e-01   -0.693  0.48823    
## factor(income)7 -7.478e-02  2.040e-01   -0.367  0.71402    
## factor(income)8 -1.467e-01  2.030e-01   -0.723  0.46988    
## factor(exint0)1 -9.559e-02  1.346e-01   -0.710  0.47768    
## factor(exint0)2 -1.855e-01  1.412e-01   -1.314  0.18914    
## factor(exint0)3 -1.542e-01  1.616e-01   -0.954  0.34007    
## block0          -4.353e-04  5.601e-04   -0.777  0.43708    
## kcal0            2.480e-05  2.964e-05    0.837  0.40285    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  2063 
## Degrees of Freedom for the fit:  23
##       Residual Deg. of Freedom:  2040 
##                       at cycle:  2 
##  
## Global Deviance:     2587.189 
##             AIC:     2633.189 
##             SBC:     2762.723 
## ******************************************************************

plot(modelo3)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.005279804 
##                        variance   =  1.02484 
##                coef. of skewness  =  0.01142201 
##                coef. of kurtosis  =  3.04374 
## Filliben correlation coefficient  =  0.9996452 
## ******************************************************************

modelo4 = update(modelo2, family = BI(mu.link = "cloglog"))
summary(modelo4)

## Warning in summary.gamlss(modelo4): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = sbp > 140 ~ initdate + factor(gender) + age +  
##     weight + weight50 + grade + factor(arth) + pkyrs + factor(diab) +  
##     factor(income) + factor(exint0) + block0 + kcal0, family = BI(mu.link = "cloglog"),  
##     data = na.omit(dados), control = gamlss.control(n.cyc = 2000,  
##         trace = FALSE)) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  cloglog
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -1.525e+02  5.687e+01  -2.681  0.00740 ** 
## initdate         9.958e-04  3.829e-04   2.601  0.00936 ** 
## factor(gender)1  6.450e-03  1.045e-01   0.062  0.95077    
## age              4.875e-02  7.603e-03   6.412 1.78e-10 ***
## weight           6.230e-03  2.220e-03   2.806  0.00507 ** 
## weight50        -6.123e-03  2.637e-03  -2.322  0.02034 *  
## grade            1.728e-04  9.383e-03   0.018  0.98531    
## factor(arth)1    1.222e-01  7.604e-02   1.606  0.10832    
## pkyrs           -1.617e-04  1.528e-03  -0.106  0.91572    
## factor(diab)2    2.266e-01  1.104e-01   2.052  0.04025 *  
## factor(diab)3    4.452e-01  1.125e-01   3.959 7.79e-05 ***
## factor(income)2 -1.262e-01  2.681e-01  -0.471  0.63783    
## factor(income)3 -5.442e-02  2.519e-01  -0.216  0.82901    
## factor(income)4  8.548e-02  2.432e-01   0.351  0.72532    
## factor(income)5  1.070e-01  2.400e-01   0.446  0.65583    
## factor(income)6 -1.165e-01  2.475e-01  -0.470  0.63806    
## factor(income)7 -3.711e-02  2.583e-01  -0.144  0.88579    
## factor(income)8 -1.439e-01  2.581e-01  -0.557  0.57731    
## factor(exint0)1 -1.262e-01  1.640e-01  -0.770  0.44148    
## factor(exint0)2 -2.418e-01  1.712e-01  -1.413  0.15783    
## factor(exint0)3 -1.874e-01  1.997e-01  -0.938  0.34832    
## block0          -4.064e-04  7.356e-04  -0.552  0.58070    
## kcal0            2.403e-05  2.431e-05   0.989  0.32301    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  2063 
## Degrees of Freedom for the fit:  23
##       Residual Deg. of Freedom:  2040 
##                       at cycle:  2 
##  
## Global Deviance:     2590.133 
##             AIC:     2636.133 
##             SBC:     2765.667 
## ******************************************************************

plot(modelo4)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.02481446 
##                        variance   =  0.9818148 
##                coef. of skewness  =  0.0609024 
##                coef. of kurtosis  =  3.000623 
## Filliben correlation coefficient  =  0.9996231 
## ******************************************************************

Tivemos umadeficiência na apresentação do resumo deste modelo, significa que não é confiável o resultado mostrado. Não consideramos então este modelo como possibilidade de ajuste.

Observando a qualidade dos resíduos decidimos escolher como mais adequado o modelo no objeto R modelo3. Percebemos ainda que muitas das covariáveis não são sgnificativamente influêntes. Procedemos então a deconsiderá-las no modelo.

modelo5 = gamlss(sbp>140 ~ initdate + age + weight + weight50 + factor(diab), 
                 data = na.omit(dados), family = BI(mu.link = "probit"),
                 control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(modelo5)

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = sbp > 140 ~ initdate + age + weight +  
##     weight50 + factor(diab), family = BI(mu.link = "probit"),  
##     data = na.omit(dados), control = gamlss.control(n.cyc = 2000,  
##         trace = FALSE)) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  probit
## Mu Coefficients:
##                 Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)   -1.128e+02  4.699e-01 -240.056  < 2e-16 ***
## initdate       7.349e-04  7.466e-06   98.428  < 2e-16 ***
## age            4.160e-02  5.830e-03    7.135 1.33e-12 ***
## weight         5.272e-03  1.756e-03    3.002  0.00272 ** 
## weight50      -4.913e-03  1.857e-03   -2.645  0.00822 ** 
## factor(diab)2  1.692e-01  8.612e-02    1.965  0.04953 *  
## factor(diab)3  3.846e-01  9.308e-02    4.132 3.73e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  2063 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  2056 
##                       at cycle:  2 
##  
## Global Deviance:     2602.208 
##             AIC:     2616.208 
##             SBC:     2655.632 
## ******************************************************************

plot(modelo5)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.02110387 
##                        variance   =  1.034507 
##                coef. of skewness  =  -0.01198588 
##                coef. of kurtosis  =  2.896262 
## Filliben correlation coefficient  =  0.9994548 
## ******************************************************************

Vejamos a significância deste modelos em relação ao modelo3.

LR.test(modelo5,modelo3)

##  Likelihood Ratio Test for nested GAMLSS models. 
##  (No check whether the models are nested is performed). 
##  
##        Null model: deviance= 2602.208 with  7 deg. of freedom 
##  Altenative model: deviance= 2587.189 with  23 deg. of freedom 
##  
##  LRT = 15.01973 with 16 deg. of freedom and p-value= 0.5231937

Significa que a resposta obtida por ambos modelos são semelhantes; sendo o modelo5 aninhado no modelo3 e, portanto, mais simples decidimos escolher este último como mais adequado.

Resultados

Lembrando que utilizamos somente os dados sem observalçoes faltantes no ajuste dos modelo, devemos agora limpar a base de dados para a apresentação dos resultados.

dados = na.omit(dados)
resposta = c(rep("Observado", length(dados$sbp)), rep("Predito", length(dados$sbp)))
resultados = data.frame(sbp = c(I(dados$sbp>140),I(fv(modelo5)>=0.5)+0.1), initdate = rep(dados$initdate,2), age = rep(dados$age,2), weight = rep(dados$weight,2), weight50 = rep(dados$weight50,2), diab = rep(dados$diab,2), resposta = resposta)

xyplot(sbp ~ initdate | diab, data = resultados, groups = resposta,
       auto.key=list(space="top", columns=2, title="Resposta", cex.title=1))

xyplot(sbp ~ age | diab, data = resultados, groups = resposta,
       auto.key=list(space="top", columns=2, title="Resposta", cex.title=1))

xyplot(sbp ~ weight | diab, data = resultados, groups = resposta,
       auto.key=list(space="top", columns=2, title="Resposta", cex.title=1))

xyplot(sbp ~ weight50 | diab, data = resultados, groups = resposta,
       auto.key=list(space="top", columns=2, title="Resposta", cex.title=1))

Vejamos a qualidade de nossas estimativas.

table(I(dados$sbp>140))

## 
## FALSE  TRUE 
##  1324   739

table(I(dados$sbp>140),I(fv(modelo5)>=0.5))

##        
##         FALSE TRUE
##   FALSE  1235   89
##   TRUE    629  110

table(I(dados$sbp>140),I(fv(modelo5)>=0.6))

##        
##         FALSE TRUE
##   FALSE  1303   21
##   TRUE    715   24

table(I(dados$sbp>140),I(fv(modelo5)>=0.4))

##        
##         FALSE TRUE
##   FALSE   997  327
##   TRUE    433  306

Temos um problema típico do modelo logístico. O quê estamos interessados em prever?