W12 in-class exercise1：Nonlinear mixed-effects modeling of MNREAD data
Ching-Fang Wu
Mon Dec 07 06:42:50 2020
1 Introduction
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
2 Data management
#set digit and stars options
options(digits = 4, show.signif.stars = FALSE)3 set path to working directory
setwd(“C:/Users/Ching-Fan Sheu/Desktop/mnread”)
# regular contrast polarity (i.e., black on white)
dta <- read.table("C:/Users/Ching-Fang Wu/Documents/data/mnread0.txt", header = T)head(dta,8)## ID pSize logRS
## 1 S101 0.0 0.9512
## 2 S101 0.1 2.0276
## 3 S101 0.2 2.0924
## 4 S101 0.3 2.2194
## 5 S101 0.4 2.2924
## 6 S101 0.5 2.2206
## 7 S101 0.6 2.2882
## 8 S101 0.7 2.2231# load packages for later use
pacman::p_load(tidyverse, nlme)# 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)")
42個受試者，依據Subject ID呈現每個人Print size與Reading time之間的關係。
# 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)")
把全部受試者Reading time曲線畫在一起，顯示出個體間的差異，顯示需使用mixed effect model。
# 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)
42個受試者，依據斜率呈現每個人Print size與Reading time之間的關係。
# try nls for all subjects
m0_g <- nlsList(SSasympOff, data = dta_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 gradient2 times caught the same error in nls(y ~ cbind(1, exp(-exp(lrc) * x)), data = xy, algorithm = “plinear”, start = list(lrc = lrc)): singular gradient
# results of parameter estimates
coef(m0_g)## Asym lrc c0
## S128 2.131 1.590 -0.30739
## S132 2.186 2.388 -0.40812
## S102 2.154 2.350 -0.04337
## S131 2.193 1.551 -0.60031
## S122 2.166 2.270 -0.34818
## S109 2.202 3.196 -0.01809
## S108 2.171 2.353 -0.23285
## S130 2.187 2.332 -0.31713
## S106 2.184 2.372 -0.38127
## S117 2.212 1.819 -0.22989
## S136 2.275 2.769 -0.08696
## S103 NA NA NA
## S113 2.240 2.177 -0.21904
## S101 2.251 2.758 -0.03504
## S123 2.274 1.853 -0.30819
## S133 2.225 2.362 -0.44646
## S112 2.270 2.075 -0.23090
## S142 2.258 2.169 -0.33552
## S118 2.267 2.607 -0.32027
## S104 2.236 2.431 -0.24719
## S120 2.268 1.521 -0.31138
## S105 2.289 2.420 -0.25797
## S124 2.290 1.788 -0.16797
## S129 2.258 2.130 -0.32087
## S134 2.290 2.236 -0.37175
## S135 2.287 2.071 -0.29440
## S125 2.272 1.961 -0.28334
## S126 2.317 1.685 -0.47288
## S141 2.323 2.082 -0.29600
## S111 2.329 2.233 -0.42614
## S115 2.306 1.623 -0.56310
## S110 2.288 2.401 -0.14867
## S107 2.321 1.720 -0.37909
## S114 2.334 2.406 -0.24746
## S139 NA NA NA
## S137 2.334 2.622 -0.26045
## S140 2.352 2.190 -0.22763
## S138 2.325 2.358 -0.37132
## S127 2.316 1.982 -0.32574
## S119 2.386 1.902 -0.30892
## S116 2.412 1.710 -0.41221
## S121 2.456 2.476 -0.22895# 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 -1324 694.7
##
## Random effects:
## Formula: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
## Level: ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## Asym 0.06717 Asym lrc
## lrc 0.33175 -0.032
## c0 0.10966 -0.103 0.492
## Residual 0.06186
##
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
## Value Std.Error DF t-value p-value
## Asym 2.2700 0.01079 626 210.47 0
## lrc 2.2028 0.05630 626 39.13 0
## c0 -0.2856 0.01739 626 -16.43 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# 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)))# model comparison
anova(m1_g, m2_g, m1a_g)## Model df AIC BIC logLik Test L.Ratio p-value
## m1_g 1 10 -1369.5 -1324.4 694.7
## m2_g 2 8 -1373.0 -1337.0 694.5 1 vs 2 0.4 0.8178
## m1a_g 3 7 -551.2 -519.7 282.6 2 vs 3 823.8 <.0001m1a_g比m2_g好
m1_g和m2_g差不多
# residuals
plot(m1a_g, id = .05, cex = .6, adj=-.05)
殘差圖結果看起來fit不是很好。
# random effects
plot(ranef(m1a_g, augFrame = T))
# parameter estimates CIs
intervals(m1a_g)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## Asym 2.2820 2.3055 2.3289
## lrc 1.4610 1.5918 1.7226
## c0 -0.4176 -0.4096 -0.4016
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: ID
## lower est. upper
## sd(Asym) 0.04262 0.06013 0.08484
## sd(lrc) 0.31114 0.39712 0.50684
## cor(Asym,lrc) -0.69649 -0.42048 -0.03608
##
## Within-group standard error:
## lower est. upper
## 0.1336 0.1415 0.1498Within-group standard error: lower est. upper 0.1336 0.1415 0.1498
# plot model fits
dta_m1a <- data.frame(dta, fit = fitted(m1a_g), fit1 = fitted(m1_g))# individual profiles on separate panel
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)")
# 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)")