Reading speed example
Szu-Yu Chen
04 Dec 2020
dtaCE1 <- read.table("C:/Users/ASUS/Desktop/data/mnread0.txt", h=T)
head(dtaCE1)## 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.22060names(dtaCE1) <- c("ID","Size","Speed")
head(dtaCE1)## ID Size Speed
## 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.22060load packages
pacman::p_load(tidyverse, nlme)1 Ggplot
1.1 Each subject
Seems like most of the subjects has a lower Reading speed when print size is < 0.0
there are still subjects who are able to remain their reading speed regardless of small print size. (e.g., S115, S131)
ot <- theme_set(theme_bw())
ggplot(dtaCE1, aes(Size, Speed)) +
geom_point(size = rel(.8), pch = 20)+
geom_line()+ facet_wrap( ~ ID) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
1.2 All subjects
Subjects showed similar response patterns for print size and reading speed.
ggplot(dtaCE1, aes(Size, Speed, group = ID)) +
geom_point(pch = 20)+
geom_line(alpha = .3) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
## Nested in subjects coerce data into grouped data object
dta_g <- groupedData(Speed ~ Size | ID, data = as.data.frame(dtaCE1),
labels = list(x = "Print size", y = "Reading speed"),
units = list(x = "(logMAR)", y = "(logWPM)") )
plot(dta_g, pch = 20, aspect = .89)
2 Models
# 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 gradient# results of parameter estimates
coef(m0_g)## 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.228950092.1 m1_g, m1a_g, m2_g
# random effects model equivalent
m1_g <- nlme(m0_g)
# show summary
summary(m1_g)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Speed ~ SSasympOff(Size, Asym, lrc, c0)
## Data: dta_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: 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)))2.2 model comparison
m1a_g is comparatively desirable
anova(m1_g, m2_g, m1a_g)## 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 <.00012.3 Residual
# residuals
plot(m1a_g, id = .05, cex = .6, adj=-.05)
# 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.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.1497979# plot model fits
dta_m1a <- data.frame(dtaCE1, fit = fitted(m1a_g), fit1 = fitted(m1_g))#
ggplot(dta_m1a, aes(Size, fit)) +
geom_point(aes(Size, Speed), pch = 1)+
geom_line(col ="steelblue", lwd = 1)+
geom_line(aes(Size, 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(Size, fit, group = ID)) +
geom_point(aes(Size, Speed), pch = 1, alpha = .6)+
geom_line(alpha = .3, col="steelblue")+
geom_line(aes(Size, fit1), col = "indianred", alpha = .3) +
labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")
The end