MLGExemplo02
Fernando L Pérez
24/05/2021
Dados
Predição de altura entre adolescentes com ESRD: Neste problema iremos “procurar” um modelo para predizer altura entre adolescentes com ESRD. O conjunto de dados pedgrowth1.txt contém dados transversais de altura (cm) e idade (meses) em N = 1.200 adolescentes norte-americanos com doença renal em estágio terminal (ESRD).
As medições foram todas feitas em 1º de janeiro, 1990 (mais ou menos um ou dois dias). Além da altura e idade, as seguintes covariáveis estão disponíveis:
• state: declara o estado no qual o paciente morava no momento da medição • female: 1 se o paciente for feminino e 0 se masculino • racegrp: Raça do paciente; 1 = Caucasiano, 2 = Afro-americano, 3 = Outro • duresrd: Quantidade de tempo que o paciente teve ESRD antes da medição (em meses) • gn: 1 se a glomerulonefrite foi a causa primária da insuficiência renal e 0 caso contrário
dados = read.table("http://www.ics.uci.edu/~dgillen/STAT211/Data/pedgrowth1.txt", header = T)
head(dados)## id state age female racegrp duresrd gn height
## 1 316235 Texas 6 1 1 0 0 59.900
## 2 404686 Minnesota 205 0 1 100 0 134.800
## 3 516093 Florida 203 1 2 0 1 162.560
## 4 21150 Indiana 185 1 1 0 0 160.020
## 5 414400 Tennessee 178 0 1 32 0 156.200
## 6 280299 Michigan 9 0 1 6 0 70.358Construa um modelo de regressão para prever a altura (height).
Estudo descritivo
library(lattice)
xyplot(height ~ age | racegrp, groups = gn, data = dados,
auto.key=list(space="top", columns=2, title="Glomerulonefrite foi a causa primária?", cex.title=1))
Observamos acima as alturas (height) segundo a idade (age) em meses nos três grupos raciais considerados. A esquerda o grupo de adolescentes Caucasianos, no centro Afro-americanos e a esquerda Outros.
min(dados$age); max(dados$age); max(dados$age)/12## [1] 0## [1] 237## [1] 19.75min(dados$height); max(dados$height)## [1] 50.8## [1] 195.58Temos um conjunto de observações em adolescentes entre 0 e quase 20 anos. Também notamos que as alturas variam entre 50cm e quase 2 metros de altura.
Perguntamos agora, o tempo que o adolescente teve a doença antes da medição influencia na resposta?
xyplot(height ~ duresrd | female, groups = gn, data = dados,
auto.key=list(space="top", columns=2, title="Glomerulonefrite foi a causa primária?", cex.title=1))
Enquanto percebemos que a idade do adolescente influencia na resposta nem o sexo, nem o tempo que o adolescente teve a doença antes da medição (duresrd) parece influenciar claramente.
Modelos
Considerando que as alturas não devem satisfazer um comportamento normal, assumimos agora que a distribuição da resposta deve ser GAMA.
library(gamlss)## Warning: package 'gamlss' was built under R version 4.0.4## Loading required package: splines## Loading required package: gamlss.data## Warning: package 'gamlss.data' was built under R version 4.0.4##
## Attaching package: 'gamlss.data'## The following object is masked from 'package:datasets':
##
## sleep## Loading required package: gamlss.dist## Warning: package 'gamlss.dist' was built under R version 4.0.4## Loading required package: MASS## Loading required package: nlme## Loading required package: parallel## ********** GAMLSS Version 5.3-2 **********## For more on GAMLSS look at https://www.gamlss.com/## Type gamlssNews() to see new features/changes/bug fixes.modeloa = gamlss(height ~ age + racegrp + gn + duresrd + female, data = dados, family = GAF(),
control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(modeloa)## ******************************************************************
## Family: c("GAF", "gamma family")
##
## Call: gamlss(formula = height ~ age + racegrp + gn + duresrd +
## female, family = GAF(), data = 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) 4.422e+00 1.336e-02 331.022 < 2e-16 ***
## age 3.302e-03 6.313e-05 52.299 < 2e-16 ***
## racegrp 7.629e-03 6.271e-03 1.217 0.22401
## gn 1.444e-02 6.855e-03 2.106 0.03543 *
## duresrd -9.631e-04 9.058e-05 -10.632 < 2e-16 ***
## female -1.987e-02 6.266e-03 -3.172 0.00156 **
## ---
## 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.8474 0.5547 -1.528 0.127
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4245 0.2265 6.289 4.5e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 1179
## Degrees of Freedom for the fit: 8
## Residual Deg. of Freedom: 1171
## at cycle: 1174
##
## Global Deviance: 9566.2
## AIC: 9582.2
## SBC: 9622.779
## ******************************************************************resposta = c(rep("Observado", length(dados$height)), rep("Predito", length(dados$height)))
resultados = data.frame(height = c(dados$height,fitted(modeloa)), age = rep(dados$age,2), racegrp = rep(dados$racegrp,2), gn = rep(dados$gn,2), duresrd = rep(dados$duresrd,2), female = rep(dados$female,2), resposta = resposta)
xyplot(height ~ age | racegrp, groups = resposta, data = resultados,
auto.key=list(space="top", columns=2, title="", cex.title=1))
Observamos que as alturas em idades baixas são completamente erráditas. Vejamos o que aconte se a distribuição da respoat for NORMAL.
modelob = gamlss(height ~ age + racegrp + gn + duresrd + female, data = dados, family = NO(),
control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(modelob)## ******************************************************************
## Family: c("NO", "Normal")
##
## Call: gamlss(formula = height ~ age + racegrp + gn + duresrd +
## female, family = NO(), data = 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) 77.392128 1.426778 54.243 < 2e-16 ***
## age 0.416552 0.006466 64.426 < 2e-16 ***
## racegrp 1.118107 0.754809 1.481 0.138792
## gn 1.834183 0.847100 2.165 0.030570 *
## duresrd -0.148689 0.010826 -13.734 < 2e-16 ***
## female -2.976268 0.764012 -3.896 0.000104 ***
## ---
## 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) 2.54449 0.02059 123.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 1179
## Degrees of Freedom for the fit: 7
## Residual Deg. of Freedom: 1172
## at cycle: 2
##
## Global Deviance: 9345.767
## AIC: 9359.767
## SBC: 9395.274
## ******************************************************************Novamente o grupo étnico não influencia na resposta. Vejamos as estimativas.
resultados1 = data.frame(height = c(dados$height,fitted(modelob)), age = rep(dados$age,2), racegrp = rep(dados$racegrp,2), gn = rep(dados$gn,2), duresrd = rep(dados$duresrd,2), female = rep(dados$female,2), resposta = resposta)
xyplot(height ~ age | gn, groups = resposta, data = resultados1,
auto.key=list(space="top", columns=2, title="", cex.title=1))
Igual acontece nesta situação. Um detalhe que percebemos é que, no modelo com resposta GAMA a relação entre a idade e a altura parece não linear.
modeloc = gamlss(height ~ age*duresrd + gn + female, data = dados, family = GAF(),
control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(modeloc)## ******************************************************************
## Family: c("GAF", "gamma family")
##
## Call: gamlss(formula = height ~ age * duresrd + gn + female,
## family = GAF(), data = 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) 4.330e+00 1.327e-02 326.160 < 2e-16 ***
## age 3.893e-03 7.862e-05 49.514 < 2e-16 ***
## duresrd 2.163e-03 3.118e-04 6.938 6.58e-12 ***
## gn 1.839e-02 6.733e-03 2.731 0.00641 **
## female -1.668e-02 6.053e-03 -2.756 0.00594 **
## age:duresrd -1.794e-05 1.745e-06 -10.279 < 2e-16 ***
## ---
## 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) -2.2947 0.5026 -4.565 5.51e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9997 0.2053 9.742 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 1179
## Degrees of Freedom for the fit: 8
## Residual Deg. of Freedom: 1171
## at cycle: 2
##
## Global Deviance: 9474.276
## AIC: 9490.276
## SBC: 9530.855
## ******************************************************************resultados2 = data.frame(height = c(dados$height,fitted(modeloc)), age = rep(dados$age,2), racegrp = rep(dados$racegrp,2), gn = rep(dados$gn,2), duresrd = rep(dados$duresrd,2), female = rep(dados$female,2), resposta = resposta)
xyplot(height ~ age | gn, groups = resposta, data = resultados2,
auto.key=list(space="top", columns=2, title="", cex.title=1))
Consideramos que desta maneira consideram-se melhor a dispersão dos dados conforme a idade aumenta. Umdetalhe importante: sabe-se que o crescimento fica repressado até uma dterminade idade e psteriormente aconte rapidamente. Ainda sabemos que as meninas crescem, em idades mais tempranas do que os meninso; mas o crescimento destes acontece mais célere neles do que nas meninas.
modelod = glm(height ~ age*duresrd + gn + female, data = dados, family = Gamma)
modeloe = segmented::segmented(modelod, seg.Z = ~age, psi = 50)
summary(modeloe)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.glm(obj = modelod, seg.Z = ~age, psi = 50)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.age 96.316 2.443
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.418e-02 1.745e-04 81.259 < 2e-16 ***
## age -6.430e-05 2.540e-06 -25.314 < 2e-16 ***
## duresrd 3.044e-06 2.579e-06 1.181 0.2380
## gn -1.057e-04 4.140e-05 -2.554 0.0108 *
## female 1.545e-04 3.824e-05 4.039 5.71e-05 ***
## U1.age 4.644e-05 2.730e-06 17.012 NA
## age:duresrd 2.751e-08 1.375e-08 2.000 0.0457 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.007898825)
##
## Null deviance: 55.570 on 1178 degrees of freedom
## Residual deviance: 10.539 on 1171 degrees of freedom
## AIC: 9340.1
##
## Convergence attained in 2 iter. (rel. change 1.8321e-06)resultados3 = data.frame(height = c(dados$height,fitted(modeloe)), age = rep(dados$age,2), racegrp = rep(dados$racegrp,2), gn = rep(dados$gn,2), duresrd = rep(dados$duresrd,2), female = rep(dados$female,2), resposta = resposta)
xyplot(height ~ age | gn, groups = resposta, data = resultados3,
auto.key=list(space="top", columns=2, title="", cex.title=1))
modelof = glm(height ~ age*duresrd + I(age^2) + gn + female, data = dados, family = Gamma)
modelog = segmented::segmented(modelof, seg.Z = ~age, psi = 50)
summary(modeloe)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.glm(obj = modelod, seg.Z = ~age, psi = 50)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.age 96.316 2.443
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.418e-02 1.745e-04 81.259 < 2e-16 ***
## age -6.430e-05 2.540e-06 -25.314 < 2e-16 ***
## duresrd 3.044e-06 2.579e-06 1.181 0.2380
## gn -1.057e-04 4.140e-05 -2.554 0.0108 *
## female 1.545e-04 3.824e-05 4.039 5.71e-05 ***
## U1.age 4.644e-05 2.730e-06 17.012 NA
## age:duresrd 2.751e-08 1.375e-08 2.000 0.0457 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.007898825)
##
## Null deviance: 55.570 on 1178 degrees of freedom
## Residual deviance: 10.539 on 1171 degrees of freedom
## AIC: 9340.1
##
## Convergence attained in 2 iter. (rel. change 1.8321e-06)resultados4 = data.frame(height = c(dados$height,fitted(modelog)), age = rep(dados$age,2), racegrp = rep(dados$racegrp,2), gn = rep(dados$gn,2), duresrd = rep(dados$duresrd,2), female = rep(dados$female,2), resposta = resposta)
xyplot(height ~ age | gn, groups = resposta, data = resultados4,
auto.key=list(space="top", columns=2, title="", cex.title=1))
Melhorou bastante, foi corrigido o afastamento em baixas idades. Como selecionar modelos? e os resíduos?