Inclass: Nonlinear mixed models
Chang-Chun Chen
2020-12-03
1 Inclass 1
1.1 Info
Problem: Reproduce the results of analysis reported in the reading time and character size example.
Data: 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(logWPM). Print size was measured in logMAR (related to degrees in visual angle).
Column 1: Subject ID
Column 2: Print size measured in logMAR
Column 3: Reading time in logWPM
1.2 Data management
# Cheung, S.-H., Kallie, C. S., Legge, G. E., & Cheong, A. M. Y. (2008).
# Nonlinear mixed-effects modeling of MNREAD data.
# Investigative Ophthalmology and Visual Science 49, 828-835.
# MNREAD data sets
# two data files:
# 1. regular contrast polarity - cheung_iovs_2008_revPol0.txt
# 2. reversed contrast polarity - cheung_iovs_2008_revPol1.txt
# Columns:
# 1. subID - subject ID 1-42
# 2. pSize - print size in logMAR
# 3. logRS - reading time in logWPM
# load packages for later use
pacman::p_load(tidyverse, nlme)
# input data
dta1 <- read.table("mnread0.txt", header = T)
names(dta1) <- c("ID","PrintSize","ReadingSpeed")
str(dta1)## 'data.frame': 670 obs. of 3 variables:
## $ ID : chr "S101" "S101" "S101" "S101" ...
## $ PrintSize : num 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
## $ ReadingSpeed: num 0.951 2.028 2.092 2.219 2.292 ...1.3 Curve Plot
1.3.1 Individual facet
Print size to Reading speed curves (individual facet)
# set plot theme to black-n-white
ot <- theme_set(theme_bw())
# individual profiles on separate panels
ggplot(dta1, aes(PrintSize, ReadingSpeed)) +
geom_point(size = rel(.8), pch = 20)+
geom_line()+
facet_wrap( ~ ID) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
When print size (logMAR) around 0.0~0.5, the reading speed (logWPM) is nearly approach the maximum.
1.4 groupedData
1.5 Lattice Graphs
Basically same as 1.3.1 plot, the order is different.
1.6 Nonlinear Mixed-Effects Models
1.6.1 SSasympOff
# try nls for all subjects
# SSasympOff: Self-Starting Nls Asymptotic Regression Model with an Offset
m0_g <- nlsList(SSasympOff, data = dta1_g)## Warning: 2 times caught the same error in nls(y ~ cbind(1, exp(-exp(lrc) * x)),
## data = xy, algorithm = "plinear", start = list(lrc = lrc)): singular gradient## Call:
## Model: ReadingSpeed ~ SSasympOff(PrintSize, Asym, lrc, c0) | ID
## Data: dta1_g
##
## Coefficients:
## Asym lrc c0
## S128 2.130726 1.590314 -0.30738633
## S132 2.186141 2.388438 -0.40811961
## S102 2.154035 2.349729 -0.04337340
## S131 2.193344 1.550622 -0.60030568
## S122 2.166236 2.269673 -0.34817782
## S109 2.202293 3.195539 -0.01808669
## S108 2.171135 2.353379 -0.23285328
## S130 2.187236 2.331566 -0.31712538
## S106 2.184381 2.372384 -0.38126823
## S117 2.212425 1.819404 -0.22988684
## S136 2.274693 2.768561 -0.08695735
## S103 NA NA NA
## S113 2.240347 2.176521 -0.21904219
## S101 2.251113 2.758345 -0.03503635
## S123 2.273536 1.852513 -0.30818832
## S133 2.225466 2.361778 -0.44646166
## S112 2.270355 2.074567 -0.23090480
## S142 2.258478 2.169279 -0.33551595
## S118 2.266828 2.606935 -0.32026871
## S104 2.235775 2.430587 -0.24718524
## S120 2.267755 1.520534 -0.31138404
## S105 2.288838 2.419815 -0.25797091
## S124 2.290248 1.787533 -0.16796544
## S129 2.257874 2.130162 -0.32087343
## S134 2.289745 2.235937 -0.37175303
## S135 2.286928 2.071451 -0.29439688
## S125 2.271665 1.961469 -0.28334236
## S126 2.316512 1.685303 -0.47288358
## S141 2.323487 2.082271 -0.29600493
## S111 2.328848 2.233323 -0.42614458
## S115 2.305744 1.623480 -0.56310360
## S110 2.287575 2.401334 -0.14866755
## S107 2.321140 1.720364 -0.37908521
## S114 2.333647 2.406131 -0.24746280
## S139 NA NA NA
## S137 2.333884 2.622207 -0.26045297
## S140 2.351785 2.189624 -0.22762796
## S138 2.324682 2.358429 -0.37132184
## S127 2.316275 1.982487 -0.32574300
## S119 2.386139 1.901716 -0.30891830
## S116 2.412031 1.710292 -0.41221119
## S121 2.455914 2.476338 -0.22895009
##
## Degrees of freedom: 638 total; 518 residual
## Residual standard error: 0.062621252 times caught the same error (S103, S139)
## Asym lrc c0
## S128 2.130726 1.590314 -0.30738633
## S132 2.186141 2.388438 -0.40811961
## S102 2.154035 2.349729 -0.04337340
## S131 2.193344 1.550622 -0.60030568
## S122 2.166236 2.269673 -0.34817782
## S109 2.202293 3.195539 -0.01808669
## S108 2.171135 2.353379 -0.23285328
## S130 2.187236 2.331566 -0.31712538
## S106 2.184381 2.372384 -0.38126823
## S117 2.212425 1.819404 -0.22988684
## S136 2.274693 2.768561 -0.08695735
## S103 NA NA NA
## S113 2.240347 2.176521 -0.21904219
## S101 2.251113 2.758345 -0.03503635
## S123 2.273536 1.852513 -0.30818832
## S133 2.225466 2.361778 -0.44646166
## S112 2.270355 2.074567 -0.23090480
## S142 2.258478 2.169279 -0.33551595
## S118 2.266828 2.606935 -0.32026871
## S104 2.235775 2.430587 -0.24718524
## S120 2.267755 1.520534 -0.31138404
## S105 2.288838 2.419815 -0.25797091
## S124 2.290248 1.787533 -0.16796544
## S129 2.257874 2.130162 -0.32087343
## S134 2.289745 2.235937 -0.37175303
## S135 2.286928 2.071451 -0.29439688
## S125 2.271665 1.961469 -0.28334236
## S126 2.316512 1.685303 -0.47288358
## S141 2.323487 2.082271 -0.29600493
## S111 2.328848 2.233323 -0.42614458
## S115 2.305744 1.623480 -0.56310360
## S110 2.287575 2.401334 -0.14866755
## S107 2.321140 1.720364 -0.37908521
## S114 2.333647 2.406131 -0.24746280
## S139 NA NA NA
## S137 2.333884 2.622207 -0.26045297
## S140 2.351785 2.189624 -0.22762796
## S138 2.324682 2.358429 -0.37132184
## S127 2.316275 1.982487 -0.32574300
## S119 2.386139 1.901716 -0.30891830
## S116 2.412031 1.710292 -0.41221119
## S121 2.455914 2.476338 -0.228950091.6.2 all random?
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: ReadingSpeed ~ SSasympOff(PrintSize, Asym, lrc, c0)
## Data: dta1_g
## AIC BIC logLik
## -1369.451 -1324.378 694.7256
##
## Random effects:
## Formula: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
## Level: ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## Asym 0.06717407 Asym lrc
## lrc 0.33175139 -0.032
## c0 0.10966221 -0.103 0.492
## Residual 0.06185672
##
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
## Value Std.Error DF t-value p-value
## Asym 2.2700465 0.01078579 626 210.46635 0
## lrc 2.2027888 0.05630121 626 39.12507 0
## c0 -0.2856456 0.01738738 626 -16.42833 0
## Correlation:
## Asym lrc
## lrc -0.072
## c0 -0.111 0.514
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.6287134 -0.3900675 0.0894696 0.5609092 2.6186525
##
## Number of Observations: 670
## Number of Groups: 421.6.3 fix c0
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: ReadingSpeed ~ SSasympOff(PrintSize, Asym, lrc, c0)
## Data: dta1_g
## AIC BIC logLik
## -551.2156 -519.6646 282.6078
##
## Random effects:
## Formula: list(Asym ~ 1, lrc ~ 1)
## Level: ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## Asym 0.06013025 Asym
## lrc 0.39711507 -0.42
## Residual 0.14145189
##
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
## Value Std.Error DF t-value p-value
## Asym 2.3054738 0.01197202 626 192.57176 0
## lrc 1.5918061 0.06675053 626 23.84709 0
## c0 -0.4096055 0.00409596 626 -100.00229 0
## Correlation:
## Asym lrc
## lrc -0.430
## c0 -0.103 0.266
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.6777240 -0.3357718 0.0207600 0.3950104 4.0791407
##
## Number of Observations: 670
## Number of Groups: 421.6.4 “Asym”?
Block the term with random effects involving “Asym”?
# model without random effects for the covariance terms involving Asym
# pdBlocked: Positive-Definite Block Diagonal Matrix
summary(m2_g <- update(m1_g, random = pdBlocked(list(Asym ~ 1, lrc + c0 ~ 1))))## Nonlinear mixed-effects model fit by maximum likelihood
## Model: ReadingSpeed ~ SSasympOff(PrintSize, Asym, lrc, c0)
## Data: dta1_g
## AIC BIC logLik
## -1373.049 -1336.991 694.5245
##
## Random effects:
## Composite Structure: Blocked
##
## Block 1: Asym
## Formula: Asym ~ 1 | ID
## Asym
## StdDev: 0.0671362
##
## Block 2: lrc, c0
## Formula: list(lrc ~ 1, c0 ~ 1)
## Level: ID
## Structure: General positive-definite
## StdDev Corr
## lrc 0.33152671 lrc
## c0 0.10968102 0.49
## Residual 0.06185418
##
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
## Value Std.Error DF t-value p-value
## Asym 2.270080 0.01078011 626 210.58029 0
## lrc 2.202614 0.05626483 626 39.14726 0
## c0 -0.285682 0.01739021 626 -16.42775 0
## Correlation:
## Asym lrc
## lrc -0.044
## c0 -0.014 0.512
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.63101037 -0.38733944 0.09036423 0.55051950 2.61624889
##
## Number of Observations: 670
## Number of Groups: 421.6.5 Comparison
## Model df AIC BIC logLik Test L.Ratio p-value
## m1_g 1 10 -1369.4512 -1324.3784 694.7256
## m2_g 2 8 -1373.0489 -1336.9907 694.5245 1 vs 2 0.4022 0.8178
## m1a_g 3 7 -551.2156 -519.6646 282.6078 2 vs 3 823.8333 <.0001may prefer m1a_g?
1.6.7 function’s random effects
random effects by subjects
# random effects
plot(ranef(m1a_g, augFrame = T),
interval="confidence",
level = 0.95, pch=19,
col.seg=gray(0.5), cex.lab=.1)
### confidence intervals
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## Asym 2.282016 2.3054738 2.3289313
## lrc 1.461018 1.5918061 1.7225944
## c0 -0.417631 -0.4096055 -0.4015801
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: ID
## lower est. upper
## sd(Asym) 0.04261859 0.06013025 0.08483731
## sd(lrc) 0.31114183 0.39711507 0.50684404
## cor(Asym,lrc) -0.69649436 -0.42048249 -0.03608080
##
## Within-group standard error:
## lower est. upper
## 0.1335709 0.1414519 0.14979791.6.8 m1a_g vs m1_g
# plot model fits
dta1_m1a <- data.frame(dta1, fit = fitted(m1a_g), fit1 = fitted(m1_g))
# red: m1_g, blue: m1a_g
ggplot(dta1_m1a, aes(PrintSize, fit)) +
geom_point(aes(PrintSize, ReadingSpeed), pch = 1)+
geom_line(col ="steelblue", lwd = .6)+
geom_line(aes(PrintSize, fit1), col = "indianred", lwd = .6) +
facet_wrap( ~ ID) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
# individual profiles on one panel, red: m1_g, blue: m1a_g
ggplot(dta1_m1a, aes(PrintSize, fit, group = ID)) +
geom_point(aes(PrintSize, ReadingSpeed), pch = 1, alpha = .6)+
geom_line(alpha = .3, col="steelblue")+
geom_line(aes(PrintSize, fit1), col = "indianred", alpha = .3) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
1.7 Reference
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.
2 Inclass 2
2.1 Info
Problem: Retrace the steps of analysis provided in the cochlear implant example
Data: 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).
ci{bdots}
2.2 Data management
# load packages for later use
pacman::p_load(bdots, tidyverse, nlme)
# load data
data(ci, package="bdots")
dta2 <- ci
names(dta2) <- c("SubjectType", "ID", "Time", "Fixations", "LookType")
str(dta2)## 'data.frame': 108216 obs. of 5 variables:
## $ SubjectType: Factor w/ 2 levels "CI","NH": 2 2 2 2 2 2 2 2 2 2 ...
## $ ID : int 36 36 36 36 36 36 36 36 36 36 ...
## $ Time : int 0 4 8 12 16 20 24 28 32 36 ...
## $ Fixations : num 0 0 0 0 0 0 0 0 0 0 ...
## $ LookType : Factor w/ 4 levels "Cohort","Rhyme",..: 1 1 1 1 1 1 1 1 1 1 ...## SubjectType ID 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.3 Plot
# set plot themes to black-n-white
ot <- theme_set(theme_bw())
dta2$ID <- as.factor(dta2$ID)
ggplot(dta2, aes(Time, Fixations)) +
geom_point(size=rel(.1), pch = 1, alpha = .2) +
facet_grid(SubjectType ~ LookType)
Even for both CI and NH, only the pattern of Target is differ from the other LookType?
2.3.1 set “Target”
# select target looking time course data
dta2_t <- subset(dta2, LookType == "Target")
# order data by subject , time
dta2_t <- dta2_t[order(dta2_t$ID, dta2_t$Time, decreasing = F), ]
# have a look
head(dta2_t)## SubjectType ID 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 Target
2.4 Function form
2.5 NLME
2.5.1 m0
# 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 = dta2g)## P S C
## 24 0.5570661 0.004107530 708.8991
## 81 0.7501355 0.004581853 804.4888
## 6 0.7480049 0.004294269 801.2305
## 90 0.7795104 0.004375641 780.4372
## 2 0.7769665 0.003588425 811.9258
## 41 0.8006724 0.006199286 728.1778
## 84 0.8190108 0.004794182 792.4237
## 30 0.8258956 0.005004583 819.3894
## 16 0.8408749 0.004715134 804.7894
## 10 0.8438018 0.004998549 788.7728
## 23 0.8912782 0.005303185 924.1540
## 75 0.8928513 0.005409086 805.2864
## 79 0.9073861 0.006026177 698.7974
## 48 0.8880807 0.004348550 792.9091
## 76 0.8908168 0.005846366 761.5108
## 4 0.9082364 0.006416738 736.6784
## 74 0.9019805 0.004944855 762.6900
## 72 0.8886978 0.007259340 822.2965
## 17 0.9086118 0.005659041 788.7869
## 29 0.9400890 0.006854420 777.8651
## 19 0.9224444 0.007384147 709.1389
## 82 0.9483024 0.007749707 734.0759
## 83 0.9332485 0.009380468 761.0706
## 94 0.9352946 0.008019069 735.1653
## 22 0.9473907 0.005805325 814.3229
## 86 0.9673920 0.007940323 690.4826
## 99 0.9680781 0.007186827 817.4132
## 3 0.9716229 0.007568480 677.6515
## 68 0.6758595 0.004333482 745.8946
## 64 0.8652494 0.007812997 761.6088
## 55 0.8759203 0.007716242 711.7787
## 98 0.8860959 0.005606832 744.7802
## 100 0.8968566 0.007229720 795.8416
## 51 0.9268904 0.006000142 705.5639
## 105 0.9254068 0.007172180 685.0009
## 38 0.9375122 0.008724398 696.7203
## 40 0.9477423 0.007432924 721.8294
## 65 0.9450332 0.009253157 671.1183
## 42 0.9622080 0.008466018 622.0181
## 71 0.9539147 0.008714007 658.6536
## 63 0.9633117 0.007946641 671.7330
## 95 0.9624662 0.007309274 725.3688
## 66 0.9669359 0.008558679 777.7416
## 106 0.9630186 0.005559693 743.0009
## 56 0.9674184 0.007670599 681.6990
## 60 0.9672931 0.008303483 679.9569
## 85 0.9712944 0.007032742 690.2741
## 102 0.9636017 0.008225961 722.6302
## 36 0.9795044 0.007956703 698.7604
## 61 0.9684397 0.007888731 683.6773
## 73 0.9772586 0.008402579 673.0469
## 47 0.9811978 0.010170584 704.2346
## 53 0.9810700 0.009323564 722.3633
## 58 0.9767827 0.007507716 682.51612.5.2 m1
did not converge
2.5.3 mean parameters values
# 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.271compare to previous: P = 0.87, S = 0.004, C = 757.8
2.5.4 m2
# update to let parameters dependent on SubjectType and simplify the covariance structure to a simple separate variances for each one of the parameters
# set P=0.903, S=0.007, C=741.271
m2 <- update(m1,
fixed = list(P + S + C ~ SubjectType),
random = pdDiag(P + S + C ~ 1),
groups = ~ ID,
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: dta2g
## AIC BIC logLik
## -142249 -142166.9 71134.49
##
## Random effects:
## Formula: list(P ~ 1, S ~ 1, C ~ 1)
## Level: ID
## Structure: Diagonal
## P.(Intercept) S.(Intercept) C.(Intercept) Residual
## StdDev: 0.07685583 0.0013435 45.50949 0.01706567
##
## Fixed effects: list(P + S + C ~ SubjectType)
## Value Std.Error DF t-value p-value
## P.(Intercept) 0.8698 0.014528 26995 59.86789 0.0000
## P.SubjectTypeNH 0.0682 0.020937 26995 3.25918 0.0011
## S.(Intercept) 0.0059 0.000254 26995 23.29822 0.0000
## S.SubjectTypeNH 0.0018 0.000366 26995 4.87299 0.0000
## C.(Intercept) 773.2397 8.605732 26995 89.85171 0.0000
## C.SubjectTypeNH -66.4032 12.400933 26995 -5.35469 0.0000
## Correlation:
## P.(In) P.STNH S.(In) S.STNH C.(In)
## P.SubjectTypeNH -0.694
## S.(Intercept) 0.000 0.000
## S.SubjectTypeNH 0.000 0.000 -0.694
## C.(Intercept) 0.000 0.000 0.000 0.000
## C.SubjectTypeNH 0.000 0.000 0.000 0.000 -0.694
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.2733183 -0.7499334 -0.1167442 0.3908058 4.9355479
##
## Number of Observations: 27054
## Number of Groups: 542.5.5 m3
# 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 | ID),
control = lmeControl(opt = 'optim',
niterEM = 500,
maxIter = 1e8,
msMaxIter = 1e8))## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Fixations ~ fx(Time, P, S, C)
## Data: dta2g
## AIC BIC logLik
## -244799 -244708.8 122410.5
##
## Random effects:
## Formula: list(P ~ 1, S ~ 1, C ~ 1)
## Level: ID
## Structure: Diagonal
## P.(Intercept) S.(Intercept) C.(Intercept) Residual
## StdDev: 0.07106324 0.001398717 46.42575 0.02976105
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.9962076
## Fixed effects: list(P + S + C ~ SubjectType)
## Value Std.Error DF t-value p-value
## P.(Intercept) 0.8786 0.014316 26995 61.37094 0.0000
## P.SubjectTypeNH 0.0652 0.020563 26995 3.17250 0.0015
## S.(Intercept) 0.0057 0.000271 26995 21.22912 0.0000
## S.SubjectTypeNH 0.0018 0.000391 26995 4.48854 0.0000
## C.(Intercept) 754.6382 9.312279 26995 81.03689 0.0000
## C.SubjectTypeNH -69.0195 13.204721 26995 -5.22688 0.0000
## Correlation:
## P.(In) P.STNH S.(In) S.STNH C.(In)
## P.SubjectTypeNH -0.696
## S.(Intercept) 0.007 -0.005
## S.SubjectTypeNH -0.005 0.008 -0.692
## C.(Intercept) 0.012 -0.008 0.001 0.000
## C.SubjectTypeNH -0.008 0.008 -0.001 0.001 -0.705
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.30299483 -1.01280477 -0.47437455 -0.06979383 1.58922732
##
## Number of Observations: 27054
## Number of Groups: 54
2.6 Residual plot
# residual plot by subjectType
plot(m3, resid(., type = "p") ~ fitted(.) | SubjectType, abline = 0, pch = '.')
# residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | ID*SubjectType,
abline =0, pch = '.', cex.lab = 0.5)
2.7 Model fit Plot
2.7.1 NH
# Normal group
ggplot(subset(dta2_t, SubjectType == "NH"), aes(Time, Fixations)) +
geom_point(size = rel(.5), pch = 1, col = "skyblue") +
geom_line(aes(y = yhat, group = ID)) +
coord_cartesian(ylim = c(0, 1)) +
facet_wrap(~ ID, ncol = 7) +
labs(x = "Time (ms)", y = "Proportion of fixations") +
theme(legend.position = "NONE")
2.7.2 CI
# implant group
ggplot(subset(dta2_t, SubjectType == "CI"), aes(Time, Fixations)) +
geom_point(size = rel(.5), pch = 1, col = "pink") +
geom_line(aes(y = yhat, group = ID)) +
coord_cartesian(ylim = c(0, 1)) +
facet_wrap(~ ID, ncol = 7) +
labs(x = "Time (ms)", y = "Proportion of fixations") +
theme(legend.position = "NONE")
2.8 Reference
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.
3 Inclass 3
3.1 Info
Problem: Replicate the results of analysis reported in the recovery from coma example.
Data: 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.
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.
3.2 Data management
# load package for later use - must have pacman installed first
pacman::p_load(tidyverse, nlme, car)
#input
data(Wong, package="car") ## Warning in data(Wong, package = "car"): data set 'Wong' not found# check the data set
# ?Wong
# save a copy
dta3 <- Wong
#rename
names(dta3) <- c("ID","Day","Duration","Gender","InjuryAge","MIQ","VIQ")
# check data structure
str(dta3)## 'data.frame': 331 obs. of 7 variables:
## $ ID : int 3358 3535 3547 3592 3728 3790 3807 3808 4253 4356 ...
## $ Day : num 30 16 40 13 19 13 37 31 40 31 ...
## $ Duration : int 4 17 1 10 6 3 5 7 3 7 ...
## $ Gender : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
## $ InjuryAge: num 20.7 55.3 55.9 61.7 30.1 ...
## $ MIQ : int 87 95 95 59 67 76 74 91 115 86 ...
## $ VIQ : int 89 77 116 73 73 69 77 110 110 83 ...## ID Day Duration Gender InjuryAge MIQ VIQ
## 1 3358 30 4 Male 20.67077 87 89
## 2 3535 16 17 Male 55.28816 95 77
## 3 3547 40 1 Male 55.91513 95 116
## 4 3592 13 10 Male 61.66461 59 73
## 5 3728 19 6 Male 30.12731 67 73
## 6 3790 13 3 Male 57.06229 76 69# themeset
ot <- theme_set(theme_bw())
# plot individuals/ Generalized additive models
ggplot(dta3, aes(sqrt(Day), MIQ, 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")
The performance IQ may be improved(?) if the person after recover from coma within 30 days
# examine the piq against coma duration (sqrt)
ggplot(dta3, aes(sqrt(Duration), MIQ)) +
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")
However, the longer period of coma, the more decreasing of performance IQ
3.3 Linear model
# use this to guess the initial values for parameters
summary(lm(MIQ ~ sqrt(Duration), data = dta3))##
## Call:
## lm(formula = MIQ ~ sqrt(Duration), data = dta3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.304 -10.894 -1.726 9.332 45.342
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.3041 1.2751 71.606 < 2e-16 ***
## sqrt(Duration) -1.2891 0.3372 -3.822 0.000158 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.83 on 329 degrees of freedom
## Multiple R-squared: 0.04252, Adjusted R-squared: 0.03961
## F-statistic: 14.61 on 1 and 329 DF, p-value: 0.0001583.4 NLM
Model: The recovery of IQ for patient i following days of coma t is Yi(t) = b0i + b1i × exp(-b2i sqrt(ti)) + εi,
where b0i is the asympotic level of IQ for patient i, -b1i is the amount of recovered IQ by patient i, b2i is the rate of recovery for patient i, where log(2)/b2i is the time to half-recovery. It is also assumed that
b0i = α0 + β0sqrt(durationi) + ε0i, b1i = α1 + β1sqrt(durationi) + ε1i,
where (ε0i, ε1i) are random effects following a bivariate normal distribution.
# fit the model
m1 <- nlme(MIQ ~ b0 + b1*exp(-b2*Day), data = dta3,
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)),
control = nlmeControl(msMaxIter = 1e8, opt = c("nlminb")))## Warning in nlme.formula(MIQ ~ b0 + b1 * exp(-b2 * Day), data = dta3, fixed =
## list(b0 ~ : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: MIQ ~ b0 + b1 * exp(-b2 * Day)
## Data: dta3
## AIC BIC logLik
## 2593.361 2627.58 -1287.68
##
## Random effects:
## Formula: list(b0 ~ 1, b1 ~ 1)
## Level: ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b0.(Intercept) 13.768822 b0.(I)
## b1.(Intercept) 2.591973 -0.999
## Residual 6.736621
##
## Fixed effects: list(b0 ~ 1 + sqrt(Duration), b1 ~ 1 + sqrt(Duration), b2 ~ 1)
## Value Std.Error DF t-value p-value
## b0.(Intercept) 97.09371 2.036464 127 47.67760 0.0000
## b0.sqrt(Duration) -1.24523 0.480469 127 -2.59169 0.0107
## b1.(Intercept) -11.14457 3.208109 127 -3.47388 0.0007
## b1.sqrt(Duration) -3.24860 1.076818 127 -3.01685 0.0031
## b2 0.00825 0.001652 127 4.99546 0.0000
## Correlation:
## b0.(I) b0.(D) b1.(I) b1.(D)
## 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.332416352 -0.365770468 0.008871749 0.382782626 2.303056146
##
## Number of Observations: 331
## Number of Groups: 2003.5 Fitted Plot
# fitted lines
ggplot(dta3_m1, aes(sqrt(Day), MIQ, group = ID)) +
geom_point(pch = 20, alpha = .5) +
geom_line(aes(sqrt(Day), 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),
Day = seq(0, 1000, 20))
newdata$MIQ <- predict(m1, newdata, level = 0)
plot(MIQ ~ Day, 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), MIQ[Duration == dc]), lwd = 2)
text(-25, MIQ[Duration == dc][1], labels = dc, adj = 0.9)} )}
grid()
text(-100, 95, labels="Length of\nComa", adj = 0, cex = .7)
The longer period of coma, the lower performance IQ score, and after the person recover from coma, the person need more days to improve the performance IQ score.
3.6 Reference
Wong, P.P., Monette, G., & Weiner, N.I. (2001). Mathematical models of cognitive recovery. Brain Injury, 15, 519-530.