W12 in-class exercise3：nonlinear mixed models
Ching-Fang Wu
Mon Dec 07 07:49:27 2020
1 Introduction
The data are for 200 patients who had suffered from traumatic brain injuries resulting in comas of varying duration. After awakening from their comas, patients were administered a standard IQ test, on average, twice per patient.
Source: Wong, P.P., Monette, G., & Weiner, N.I. (2001). Mathematical models of cognitive recovery. Brain Injury, 15, 519-530.
Data: Wong{car} R
Column 1: id - patient ID number. Column 2: days - number of days post coma at which IQs were measured. Column 3: duration - duration of the coma in days. Column 4: sex - a factor with levels Female and Male. Column 5: age - in years at the time of injury. Column 6: piq - performance (i.e., mathematical) IQ. Column 7: viq - verbal IQ.
2 Data management
#
options(digits = 4, show.signif.stars = FALSE)# load package for later use - must have pacman installed first
pacman::p_load(tidyverse, nlme, car)3 set path to working directory
setwd(“C:/Users/Ching-Fan Sheu/Desktop/mlm_nlm”)
# check the data set
?Wong## starting httpd help server ... done#input
data(Wong, package="car") ## Warning in data(Wong, package = "car"): data set 'Wong' not found# save a copy
dta <- Wong# check data structure
str(dta)## 'data.frame': 331 obs. of 7 variables:
## $ id : int 3358 3535 3547 3592 3728 3790 3807 3808 4253 4356 ...
## $ days : num 30 16 40 13 19 13 37 31 40 31 ...
## $ duration: int 4 17 1 10 6 3 5 7 3 7 ...
## $ sex : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
## $ age : num 20.7 55.3 55.9 61.7 30.1 ...
## $ piq : int 87 95 95 59 67 76 74 91 115 86 ...
## $ viq : int 89 77 116 73 73 69 77 110 110 83 ...# see the first 6 lines
head(dta)## id days duration sex age piq viq
## 1 3358 30 4 Male 20.67 87 89
## 2 3535 16 17 Male 55.29 95 77
## 3 3547 40 1 Male 55.92 95 116
## 4 3592 13 10 Male 61.66 59 73
## 5 3728 19 6 Male 30.13 67 73
## 6 3790 13 3 Male 57.06 76 69#
ot <- theme_set(theme_bw())# plot individuals
ggplot(dta, aes(sqrt(days), piq, group = id)) +
geom_point(pch = 20) +
geom_line() +
stat_smooth(aes(group=1), method = "gam", formula = y ~ s(x)) +
labs(x = "Days after recovering from coma (square-root)",
y = "Performance IQ score")
# examine the piq against coma duratin (sqrt)
ggplot(dta, aes(sqrt(duration), piq)) +
geom_point(pch = 20) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) +
labs(x = "Duration of coma (square-root)", y = "Performance IQ score")
# use this to guess the initial values for parameters
summary(lm(piq ~ sqrt(duration), data = dta))##
## Call:
## lm(formula = piq ~ sqrt(duration), data = dta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.30 -10.89 -1.73 9.33 45.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.304 1.275 71.61 < 2e-16
## sqrt(duration) -1.289 0.337 -3.82 0.00016
##
## Residual standard error: 14.8 on 329 degrees of freedom
## Multiple R-squared: 0.0425, Adjusted R-squared: 0.0396
## F-statistic: 14.6 on 1 and 329 DF, p-value: 0.000158# fit the model
m1 <- nlme(piq ~ b0 + b1*exp(-b2*days), data = dta,
fixed = list(b0 ~ 1 + sqrt(duration),
b1 ~ 1 + sqrt(duration),
b2 ~ 1),
random = list(id = list(b0 ~ 1, b1 ~ 1)),
start = list(fixed = c(100, -1.5, -10, 0, 0.007)))## Warning in nlme.formula(piq ~ b0 + b1 * exp(-b2 * days), data = dta, fixed =
## list(b0 ~ : Iteration 1, LME step: nlminb() did not converge (code = 1). Do
## increase 'msMaxIter'!## Warning in nlme.formula(piq ~ b0 + b1 * exp(-b2 * days), data = dta, fixed =
## list(b0 ~ : Iteration 2, LME step: nlminb() did not converge (code = 1). Do
## increase 'msMaxIter'!Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase ‘msMaxIter’!Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase ‘msMaxIter’!
# show output summary
summary(m1)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: piq ~ b0 + b1 * exp(-b2 * days)
## Data: dta
## AIC BIC logLik
## 2593 2628 -1288
##
## Random effects:
## Formula: list(b0 ~ 1, b1 ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b0.(Intercept) 13.769 b0.(I)
## b1.(Intercept) 2.600 -0.996
## Residual 6.736
##
## Fixed effects: list(b0 ~ 1 + sqrt(duration), b1 ~ 1 + sqrt(duration), b2 ~ 1)
## Value Std.Error DF t-value p-value
## b0.(Intercept) 97.09 2.037 127 47.68 0.0000
## b0.sqrt(duration) -1.25 0.480 127 -2.59 0.0107
## b1.(Intercept) -11.15 3.208 127 -3.47 0.0007
## b1.sqrt(duration) -3.25 1.077 127 -3.02 0.0031
## b2 0.01 0.002 127 5.00 0.0000
## Correlation:
## b0.(I) b0.s() b1.(I) b1.s()
## b0.sqrt(duration) -0.724
## b1.(Intercept) -0.596 0.463
## b1.sqrt(duration) 0.463 -0.455 -0.789
## b2 -0.309 0.013 0.092 -0.380
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.332390 -0.365703 0.008973 0.382757 2.303088
##
## Number of Observations: 331
## Number of Groups: 200# augment data with fitted values
dta_m1 <- data.frame(dta, fit = predict(m1))# fitted lines
ggplot(dta_m1, aes(sqrt(days), piq, group = id)) +
geom_point(pch = 20, alpha = .5) +
geom_line(aes(sqrt(days), fit, group = id)) +
labs(x = "Days after recovering from coma (square-root)",
y = "Performance IQ score")
##
# plot mode fits - basic graphics - detail contol
#
newdata <- expand.grid(duration = c(1, 10, 20, 50, 100, 200),
days = seq(0, 1000, 20))
newdata$piq <- predict(m1, newdata, level = 0)plot(piq ~ days, type="n", data = newdata,
xlab = "Days Post Coma", ylab = "Average Performance IQ",
ylim = c(15, 105), xlim = c(-100, 1000), axes = FALSE, frame=TRUE)
axis(2)
axis(4)
axis(1, at = seq(0, 1000, by=100))
for (dc in c(1, 10, 20, 50, 100, 200)){
with(newdata, {
lines(spline(seq(0, 1000, 20), piq[duration == dc]), lwd = 2)
text(-25, piq[duration == dc][1], labels = dc, adj = 0.9)
}
)}
grid()
text(-100, 95, labels="Length of\nComa", adj = 0, cex = .7)
昏迷時間越長，智商表現越低，且清醒後需要更多時間來改善智商的表現。 # ###