Nonlinear mixed models: Exercises
Hao-Lun Fu
2020-12-07
1 Exercises 1:
Reproduce the results of analysis reported in the reading time and character size example.
MNREAD Acuity Charts (Optelec, US Inc.) were used to measure reading speed as a function of print size. MNREAD data were collected from 42 adult observers with normal vision. Reading speed was measured in log word per minute. Print size was measured in in logMAR (related to degrees in visual angle).
Source: Cheung, C., Kallie, S., Legge, G.E., & Cheong, A.M.Y. (2008). Nonlinear Mixed-Effects Modeling of MNREAD Data. Investigative Ophthalmology & Visual Science, 49(2), 828-835.
Column 1: Subject ID Column 2: Reading time in logWPM Column 3: Print size measured in logMAR
1.1 Load data file
# load packages for later use
pacman::p_load(tidyverse, nlme)
# regular contrast polarity (i.e., black on white)
dta <- read.table("Data/mnread0.txt", header = T)
head(dta) ID pSize logRS
1 S101 0.0 0.95121
2 S101 0.1 2.02760
3 S101 0.2 2.09240
4 S101 0.3 2.21940
5 S101 0.4 2.29240
6 S101 0.5 2.220601.2 Data Visualization
# set plot theme to black-n-white
ot <- theme_set(theme_bw())
# individual profiles on separate panels
ggplot(dta, aes(pSize, logRS)) +
geom_point(size = rel(.8), pch = 20)+
geom_line()+
facet_wrap( ~ ID) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
1.3 Individual profiles on one panel
ggplot(dta, aes(pSize, logRS, group = ID)) +
geom_point(pch = 20)+
geom_line(alpha = .3) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
1.4 Coerce data into groupedData object
dta_g <- groupedData(logRS ~ pSize | ID, data = as.data.frame(dta),
labels = list(x = "Print size", y = "Reading speed"),
units = list(x = "(logMAR)", y = "(logWPM)") )
# lattice plot by subject
plot(dta_g, pch = 20, aspect = .89)
1.5 NLS for all subjects
m0_g <- nlsList(SSasympOff, data = dta_g)
# results of parameter estimates
coef(m0_g) Asym lrc c0
S128 2.1307 1.5903 -0.307386
S132 2.1861 2.3884 -0.408120
S102 2.1540 2.3497 -0.043373
S131 2.1933 1.5506 -0.600306
S122 2.1662 2.2697 -0.348178
S109 2.2023 3.1955 -0.018087
S108 2.1711 2.3534 -0.232853
S130 2.1872 2.3316 -0.317125
S106 2.1844 2.3724 -0.381268
S117 2.2124 1.8194 -0.229887
S136 2.2747 2.7686 -0.086957
S103 NA NA NA
S113 2.2403 2.1765 -0.219042
S101 2.2511 2.7583 -0.035036
S123 2.2735 1.8525 -0.308188
S133 2.2255 2.3618 -0.446462
S112 2.2704 2.0746 -0.230905
S142 2.2585 2.1693 -0.335516
S118 2.2668 2.6069 -0.320269
S104 2.2358 2.4306 -0.247185
S120 2.2678 1.5205 -0.311384
S105 2.2888 2.4198 -0.257971
S124 2.2902 1.7875 -0.167965
S129 2.2579 2.1302 -0.320873
S134 2.2897 2.2359 -0.371753
S135 2.2869 2.0715 -0.294397
S125 2.2717 1.9615 -0.283342
S126 2.3165 1.6853 -0.472884
S141 2.3235 2.0823 -0.296005
S111 2.3288 2.2333 -0.426145
S115 2.3057 1.6235 -0.563104
S110 2.2876 2.4013 -0.148668
S107 2.3211 1.7204 -0.379085
S114 2.3336 2.4061 -0.247463
S139 NA NA NA
S137 2.3339 2.6222 -0.260453
S140 2.3518 2.1896 -0.227628
S138 2.3247 2.3584 -0.371322
S127 2.3163 1.9825 -0.325743
S119 2.3861 1.9017 -0.308918
S116 2.4120 1.7103 -0.412211
S121 2.4559 2.4763 -0.2289501.6 Random effects model equivalent
m1_g <- nlme(m0_g)
# show summary
summary(m1_g)Nonlinear mixed-effects model fit by maximum likelihood
Model: logRS ~ SSasympOff(pSize, Asym, lrc, c0)
Data: dta_g
AIC BIC logLik
-1369.5 -1324.4 694.73
Random effects:
Formula: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
Level: ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
Asym 0.067174 Asym lrc
lrc 0.331751 -0.032
c0 0.109662 -0.103 0.492
Residual 0.061857
Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
Value Std.Error DF t-value p-value
Asym 2.27005 0.010786 626 210.466 0
lrc 2.20279 0.056301 626 39.125 0
c0 -0.28565 0.017387 626 -16.428 0
Correlation:
Asym lrc
lrc -0.072
c0 -0.111 0.514
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-6.62871 -0.39007 0.08947 0.56091 2.61865
Number of Observations: 670
Number of Groups: 42 1.7 Model comparison
# Fix c0
m1a_g <- update(m1_g, random = Asym + lrc ~ 1)
# Model without random effects for the covariance terms involving Asym
m2_g <- update(m1_g, random = pdBlocked(list(Asym ~ 1, lrc + c0 ~ 1)))
#
anova(m1_g, m2_g, m1a_g) Model df AIC BIC logLik Test L.Ratio p-value
m1_g 1 10 -1369.45 -1324.38 694.73
m2_g 2 8 -1373.05 -1336.99 694.52 1 vs 2 0.40 0.8178
m1a_g 3 7 -551.22 -519.66 282.61 2 vs 3 823.83 <.00011.8 Plot residuals
plot(m1a_g, id = .05, cex = .6, adj=-.05)
1.9 Plot random effects
plot(ranef(m1a_g, augFrame = T))
1.10 Parameter estimates CIs
intervals(m1a_g)Approximate 95% confidence intervals
Fixed effects:
lower est. upper
Asym 2.28202 2.30547 2.32893
lrc 1.46102 1.59181 1.72259
c0 -0.41763 -0.40961 -0.40158
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: ID
lower est. upper
sd(Asym) 0.042617 0.06013 0.084841
sd(lrc) 0.311089 0.39712 0.506930
cor(Asym,lrc) -0.697018 -0.42048 -0.035064
Within-group standard error:
lower est. upper
0.13357 0.14145 0.14980 1.11 Plot model fits
dta_m1a <- data.frame(dta, fit = fitted(m1a_g), fit1 = fitted(m1a_g))
#
ggplot(dta_m1a, aes(pSize, fit)) +
geom_point(aes(pSize, logRS), pch = 1)+
geom_line(col ="steelblue", lwd = 1)+
geom_line(aes(pSize, fit1), col = "indianred", lwd = 1) +
facet_wrap( ~ ID) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
1.12 Individual profiles on one panel
ggplot(dta_m1a, aes(pSize, fit, group = ID)) +
geom_point(aes(pSize, logRS), pch = 1, alpha = .6)+
geom_line(alpha = .3, col="steelblue")+
geom_line(aes(pSize, fit1), col = "indianred", alpha = .3) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
The results revealed that the reading time increased as the print size of the word increased.
2 Exercise 2:
Retrace the steps of analysis provided in the cochlear implant example
Twenty-eight postlingually-deafened cochlear implant (CI) users and 26 normal hearing (NH) adults perform a spoken-word recognition task. On each trial, 4 pictures appear on a monitor. The subject hears the name of one of them (e.g., wizard as the target) and clicks on the referent. The competitors, in this case, are: lizard (rhyme), whistle (cohort), baggage (unrelated).
Source: Farris-Trimble, A., McMurray, B., Cigrand, N., & Tomblin, J.B. (2014). The process of spoken word recognition in the face of signal degradation: Cochlear implant users and normal-hearing listeners. Journal of Experimental Psychology: Human Perception and Performance, 40(1), 308-327.
# load packages for later use
pacman::p_load(bdots, tidyverse, nlme)
# load data
data(ci)
#
head(ci) protocol Subject Time Fixations LookType
1 NH 36 0 0 Cohort
2 NH 36 4 0 Cohort
3 NH 36 8 0 Cohort
4 NH 36 12 0 Cohort
5 NH 36 16 0 Cohort
6 NH 36 20 0 Cohort2.1 Data Visualization
# set plot themes to black-n-white
ot <- theme_set(theme_bw())
#ci$Subject <- as.factor(ci$Subject)
ggplot(ci, aes(Time, Fixations)) +
geom_point(size=rel(.1), pch = 1, alpha = .2) +
facet_grid(protocol ~ LookType)
2.2 Data manipulation
# select target looking time course data
dta_t <- subset(ci, LookType == "Target")
# order data by subject , time
dta_t <- dta_t[order(dta_t$Subject, dta_t$Time, decreasing = F), ]
# have a look
head(dta_t) protocol Subject Time Fixations LookType
39079 CI 2 0 0 Target
39080 CI 2 4 0 Target
39081 CI 2 8 0 Target
39082 CI 2 12 0 Target
39083 CI 2 16 0 Target
39084 CI 2 20 0 Target2.3 Convert to groupedData object for nlme processing
dtag <- groupedData(Fixations ~ Time | Subject, outer = ~ protocol,
data = as.data.frame(dta_t),
labels = list(x = "Time", y = "Fixations"),
units = list(x = "(ms)", y = "(proportion)") )
# lattice plot by subject
plot(dtag, pch = 20, aspect = .89)
2.4 Plot by subject given group = protocol
p <- plot(dtag, outer = TRUE, cex = .1, layout = c(2,1), aspect = .89)
update(p, legend = NULL)
2.5 Model 0
# define the 3-parameter curve fitting function
fx <- function(x, P, S, C) {
P / (1 + exp(S*(C - x)/P))
}
#
# fit nls separately for all subjects
#
m0 <- nlsList(Fixations ~ fx(Time, P, S, C),
start = c(P = 0.87, S = 0.004, C = 757.8),
data = dtag)
# results of parameter estimates
coef(m0) P S C
24 0.55707 0.0041075 708.90
81 0.75014 0.0045819 804.49
6 0.74800 0.0042943 801.23
90 0.77951 0.0043756 780.44
2 0.77697 0.0035884 811.93
41 0.80067 0.0061993 728.18
84 0.81901 0.0047942 792.42
30 0.82590 0.0050046 819.39
16 0.84087 0.0047151 804.79
10 0.84380 0.0049985 788.77
23 0.89128 0.0053032 924.15
75 0.89285 0.0054091 805.29
79 0.90739 0.0060262 698.80
48 0.88808 0.0043485 792.91
76 0.89082 0.0058464 761.51
4 0.90824 0.0064167 736.68
74 0.90198 0.0049449 762.69
72 0.88870 0.0072593 822.30
17 0.90861 0.0056590 788.79
29 0.94009 0.0068544 777.87
19 0.92244 0.0073841 709.14
82 0.94830 0.0077497 734.08
83 0.93325 0.0093805 761.07
94 0.93529 0.0080191 735.17
22 0.94739 0.0058053 814.32
86 0.96739 0.0079403 690.48
99 0.96808 0.0071868 817.41
3 0.97162 0.0075685 677.65
68 0.67586 0.0043335 745.89
64 0.86525 0.0078130 761.61
55 0.87592 0.0077162 711.78
98 0.88610 0.0056068 744.78
100 0.89686 0.0072297 795.84
51 0.92689 0.0060001 705.56
105 0.92541 0.0071722 685.00
38 0.93751 0.0087244 696.72
40 0.94774 0.0074329 721.83
65 0.94503 0.0092532 671.12
42 0.96221 0.0084660 622.02
71 0.95391 0.0087140 658.65
63 0.96331 0.0079466 671.73
95 0.96247 0.0073093 725.37
66 0.96694 0.0085587 777.74
106 0.96302 0.0055597 743.00
56 0.96742 0.0076706 681.70
60 0.96729 0.0083035 679.96
85 0.97129 0.0070327 690.27
102 0.96360 0.0082260 722.63
36 0.97950 0.0079567 698.76
61 0.96844 0.0078887 683.68
73 0.97726 0.0084026 673.05
47 0.98120 0.0101706 704.23
53 0.98107 0.0093236 722.36
58 0.97678 0.0075077 682.522.6 Model 1: Covert fitting separately by subject to subjects as random effects
m1 <- nlme(m0)
# extract mean parameters values to use as initial values
prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)
prm0 P S C
0.903 0.007 741.271 2.7 Model 2
# update to let parameters dependent on protocol
# and simplify the covariance structure to a simple separate
# variances for each one of the parameters
#
m2 <- update(m1,
fixed = list(P + S + C ~ protocol),
random = pdDiag(P + S + C ~ 1),
groups = ~ Subject,
start = list(fixed = c(P = c(prm0[1], 0),
S = c(prm0[2], 0),
C = c(prm0[3], 0))))
summary(m2)Nonlinear mixed-effects model fit by maximum likelihood
Model: Fixations ~ fx(Time, P, S, C)
Data: dtag
AIC BIC logLik
-142249 -142167 71134
Random effects:
Formula: list(P ~ 1, S ~ 1, C ~ 1)
Level: Subject
Structure: Diagonal
P.(Intercept) S.(Intercept) C.(Intercept) Residual
StdDev: 0.076856 0.0013435 45.509 0.017066
Fixed effects: list(P + S + C ~ protocol)
Value Std.Error DF t-value p-value
P.(Intercept) 0.87 0.0145 26995 59.868 0.0000
P.protocolNH 0.07 0.0209 26995 3.259 0.0011
S.(Intercept) 0.01 0.0003 26995 23.298 0.0000
S.protocolNH 0.00 0.0004 26995 4.873 0.0000
C.(Intercept) 773.24 8.6057 26995 89.852 0.0000
C.protocolNH -66.40 12.4009 26995 -5.355 0.0000
Correlation:
P.(In) P.prNH S.(In) S.prNH C.(In)
P.protocolNH -0.694
S.(Intercept) 0.000 0.000
S.protocolNH 0.000 0.000 -0.694
C.(Intercept) 0.000 0.000 0.000 0.000
C.protocolNH 0.000 0.000 0.000 0.000 -0.694
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.27332 -0.74993 -0.11674 0.39081 4.93555
Number of Observations: 27054
Number of Groups: 54 2.8 Model 3
# update model 2 to include auto-regressive dependence over time
# results have convergence issues - takes a few minutes
# set lmeControl options - change optimizer, etc.
m3 <- update(m2, corr = corAR1(0.95, form = ~ 1 | Subject),
control = lmeControl(opt = 'optim',
niterEM = 500,
maxIter = 1e6,
msMaxIter = 1e3))
# show summary
summary(m3)Nonlinear mixed-effects model fit by maximum likelihood
Model: Fixations ~ fx(Time, P, S, C)
Data: dtag
AIC BIC logLik
-244799 -244709 122411
Random effects:
Formula: list(P ~ 1, S ~ 1, C ~ 1)
Level: Subject
Structure: Diagonal
P.(Intercept) S.(Intercept) C.(Intercept) Residual
StdDev: 0.071063 0.0013987 46.426 0.029761
Correlation Structure: AR(1)
Formula: ~1 | Subject
Parameter estimate(s):
Phi
0.99621
Fixed effects: list(P + S + C ~ protocol)
Value Std.Error DF t-value p-value
P.(Intercept) 0.88 0.0143 26995 61.371 0.0000
P.protocolNH 0.07 0.0206 26995 3.173 0.0015
S.(Intercept) 0.01 0.0003 26995 21.229 0.0000
S.protocolNH 0.00 0.0004 26995 4.489 0.0000
C.(Intercept) 754.64 9.3123 26995 81.037 0.0000
C.protocolNH -69.02 13.2047 26995 -5.227 0.0000
Correlation:
P.(In) P.prNH S.(In) S.prNH C.(In)
P.protocolNH -0.696
S.(Intercept) 0.007 -0.005
S.protocolNH -0.005 0.008 -0.692
C.(Intercept) 0.012 -0.008 0.001 0.000
C.protocolNH -0.008 0.008 -0.001 0.001 -0.705
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.302904 -1.012773 -0.474290 -0.069741 1.589378
Number of Observations: 27054
Number of Groups: 54 2.9 Check the auto-correlation function of residuals
plot(ACF(m3, resType = "n"), alpha = 0.05)
2.10 Residual plot to identify outlier
plot(m3, id = 0.05, pch = '.', adj = -0.1, cex = .4)
2.11 Residual plot by protocol
plot(m3, resid(., type = "p") ~ fitted(.) | protocol, abline = 0, pch = '.')
2.12 Residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | Subject, abline = 0, pch = '.')
2.13 Residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | Subject*protocol,
abline = 0, pch = '.')
2.14 Plot fitted lines over observations by subject
plot(augPred(m3, level = 0:1))
2.15 Model Prediction
# fortify data with fitted values
dta_t$yhat <- fitted(m3)
# Normal group
ggplot(subset(dta_t, protocol == "NH"), aes(Time, Fixations)) +
geom_point(size = rel(.5), pch = 1, col = "skyblue") +
geom_line(aes(y = yhat, group = Subject)) +
coord_cartesian(ylim = c(0, 1)) +
facet_wrap(~ Subject, ncol = 7) +
labs(x = "Time (ms)", y = "Proportion of fixations") +
theme(legend.position = "NONE")
# implant group
ggplot(subset(dta_t, protocol == "CI"), aes(Time, Fixations)) +
geom_point(size = rel(.5), pch = 1, col = "pink") +
geom_line(aes(y = yhat, group = Subject)) +
coord_cartesian(ylim = c(0, 1)) +
facet_wrap(~ Subject, ncol = 7) +
labs(x = "Time (ms)", y = "Proportion of fixations") +
theme(legend.position = "NONE")
The results revealed that the protocol (cochlear implant (CI) users vs. normal hearing (NH) adults) and individual differences had a significant effect on the spoken-word recognition task, the NH adults had better sensitivity curve function than the CI users.
3 Exercise 3:
Replicate the results of analysis reported in the recovery from coma example.
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.
3.1 Load data file
# load package for later use - must have pacman installed first
pacman::p_load(tidyverse, nlme, car)
# 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.671 87 89
2 3535 16 17 Male 55.288 95 77
3 3547 40 1 Male 55.915 95 116
4 3592 13 10 Male 61.665 59 73
5 3728 19 6 Male 30.127 67 73
6 3790 13 3 Male 57.062 76 693.2 Data Visualization
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")
3.3 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")
3.4 LM for duration and IQ
# 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.0001583.5 Model 1
# 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)))
# 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.4 2627.6 -1287.7
Random effects:
Formula: list(b0 ~ 1, b1 ~ 1)
Level: id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
b0.(Intercept) 13.7691 b0.(I)
b1.(Intercept) 2.6010 -0.996
Residual 6.7362
Fixed effects: list(b0 ~ 1 + sqrt(duration), b1 ~ 1 + sqrt(duration), b2 ~ 1)
Value Std.Error DF t-value p-value
b0.(Intercept) 97.095 2.0365 127 47.676 0.0000
b0.sqrt(duration) -1.245 0.4805 127 -2.592 0.0107
b1.(Intercept) -11.145 3.2081 127 -3.474 0.0007
b1.sqrt(duration) -3.248 1.0768 127 -3.017 0.0031
b2 0.008 0.0017 127 4.996 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.3323903 -0.3656999 0.0089778 0.3827555 2.3030905
Number of Observations: 331
Number of Groups: 200 3.6 Model Prediction
# 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)
The results revealed that the duration of coma and the days after recovering from coma significantly affected the IQ performance for these patients who had suffered from traumatic brain injuries resulting in comas.