During a 3-week period, subjects played a chess board game called, “Fox & Geese”. Each subject played up to 27 games. The number of moves made by the subject before commiting a catastrophic error (resulting in defeat) was the outcome of interest. The reading scores of these subjects were also collected as a potential covariate.
Source: Tivnan, T. (1980). Dissertation at Harvard Graduate School of Education. Reported in Singer, J.D., & Willet, J.B. (2003). Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence.
Column 1: Subject ID Column 2: Game played Column 3: Number of moves made before defeat Column 4: Reading score
# load packages
pacman::p_load(tidyverse, nlme, car,lattice,nlstools)
# input data
dta<- read.table("C:/Users/Ching-Fang Wu/Documents/data//fox_geese.txt", h=T)
# show first 8 lines
knitr::kable(head(dta,8))
ID | Game | Moves | Read |
---|---|---|---|
S101 | 1 | 4 | 1.4 |
S101 | 2 | 7 | 1.4 |
S101 | 3 | 8 | 1.4 |
S101 | 4 | 3 | 1.4 |
S101 | 5 | 3 | 1.4 |
S101 | 6 | 3 | 1.4 |
S101 | 7 | 7 | 1.4 |
S101 | 8 | 6 | 1.4 |
str(dta)
## 'data.frame': 445 obs. of 4 variables:
## $ ID : chr "S101" "S101" "S101" "S101" ...
## $ Game : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Moves: int 4 7 8 3 3 3 7 6 3 7 ...
## $ Read : num 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 ...
dta$ID <- as.factor(dta$ID)
#個別受試者移動次數與遊戲之間的關係(依照Subject ID順序)
xyplot(Moves ~ Game | ID , data=dta,
type=c("g", "p"), cex=.5,
xlab="Number of game played",
ylab="Number of moves made before defeat")
多數受試者玩過的遊戲種類越多,移動次數越多,但是S101、S103和S106例外,他們的移動次數沒有太大變化。
##個別受試者移動次數與遊戲之間的關係(依照斜率高低)
dta_g <- groupedData(Moves ~ Game | ID, outer = ~ Read,
data = as.data.frame(dta),
labels = list(x = "Game", y = "Moves"),
units = list(x = "", y = "") )
plot(dta_g, pch = 1, aspect = .9)
#將全部受試者移動次數與遊戲之間的關係圖畫在一張
ggplot(dta, aes(Game, Moves, group=ID))+
geom_point(alpha=.3)+
geom_line(alpha=.5)+
labs(x="Game", "Moves")
這張圖顯示出個體間的差異。
Model: Let Yij be the number of moves for a subject i at jth game played. Yij = 1 + 19 / (1 + b0i exp(-b1iTimeij)) + εij,
where Timeij is the number of times game played by subject i. The parameters (b0i, b1i) are random effects.
b0i (150, 15, or 1.5) b1i (0.1, 0.3, or 0.5)
mfx <- function(x, A, B) { (1 + (19 /(1 + (A*exp(-B*x)))))}
m0 <- nlsList(Moves ~ mfx(Game, A, B),
start = c(A = 150, B = 0.1),
data = dta_g)
## Warning: 1 error caught in nls(model, data = data, control = controlvals, start
## = start): number of iterations exceeded maximum of 50
coef(m0)
## A B
## S113 3.444376 0.019393820
## S110 19.021257 0.083381612
## S101 3.085531 -0.034697912
## S102 28.499570 0.148899148
## S103 5.385204 0.005747663
## S104 6.257776 0.133481440
## S108 26.681234 0.250404162
## S116 7.093699 0.141653735
## S106 8.742240 0.017753536
## S105 16.633833 0.125658158
## S117 NA NA
## S115 16.674963 0.160175972
## S107 3.146240 0.047900002
## S114 101.147125 0.298797401
## S112 88.215308 0.176818526
## S111 56.277746 0.198327030
## S109 8.227554 0.155514017
##m1
m1 <- nlme(m0, random=pdDiag(A + B ~ 1))
prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)
prm0
## A B
## 10.357 0.109
##m2
m2 <- update(m1, weights=varIdent(form = ~ 1 | ID),
control=nlmeControl(maxIter=1e8, msMaxIter=1e8, niterEM=500))
summary(m2)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Moves ~ mfx(Game, A, B)
## Data: dta_g
## AIC BIC logLik
## 2456.174 2542.234 -1207.087
##
## Random effects:
## Formula: list(A ~ 1, B ~ 1)
## Level: ID
## Structure: Diagonal
## A B Residual
## StdDev: 3.468549 0.05936851 3.982061
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | ID
## Parameter estimates:
## S113 S110 S101 S102 S103 S104 S108 S116
## 1.0000000 0.8277653 0.4279184 1.1215850 0.5515780 1.1217129 0.7118691 1.2247194
## S106 S105 S117 S115 S107 S114 S112 S111
## 0.3392785 0.8255157 1.0328394 0.9281161 0.6649241 0.9527503 0.9995639 1.1471650
## S109
## 1.2466181
## Fixed effects: list(A ~ 1, B ~ 1)
## Value Std.Error DF t-value p-value
## A 8.579734 1.2977639 427 6.611167 0
## B 0.096326 0.0159375 427 6.044003 0
## Correlation:
## A
## B 0.287
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.1980973 -0.6733032 -0.1749343 0.5312265 4.0182717
##
## Number of Observations: 445
## Number of Groups: 17
##m3
m3 <- update(m2,
fixed = list(A + B ~ Read),
random = pdDiag(A + B ~ 1),
groups = ~ ID,
start = list(fixed = c(A = c(prm0[1], 0),
B = c(prm0[2], 0))))
summary(m3)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Moves ~ mfx(Game, A, B)
## Data: dta_g
## AIC BIC logLik
## 2452.096 2546.352 -1203.048
##
## Random effects:
## Formula: list(A ~ 1, B ~ 1)
## Level: ID
## Structure: Diagonal
## A.(Intercept) B.(Intercept) Residual
## StdDev: 2.026809 0.05351836 3.894973
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | ID
## Parameter estimates:
## S113 S110 S101 S102 S103 S104 S108 S116
## 1.0000000 0.8633940 0.4407039 1.1670976 0.5603937 1.1435745 0.7424576 1.2491533
## S106 S105 S117 S115 S107 S114 S112 S111
## 0.3458037 0.8407405 1.0628409 0.9404110 0.8148864 0.9537729 0.9889874 1.1371776
## S109
## 1.3115064
## Fixed effects: list(A + B ~ Read)
## Value Std.Error DF t-value p-value
## A.(Intercept) -4.631553 4.837131 425 -0.9575000 0.3389
## A.Read 8.021219 3.152020 425 2.5447863 0.0113
## B.(Intercept) -0.001241 0.040061 425 -0.0309764 0.9753
## B.Read 0.054672 0.019629 425 2.7852689 0.0056
## Correlation:
## A.(In) A.Read B.(In)
## A.Read -0.975
## B.(Intercept) 0.306 -0.273
## B.Read -0.395 0.389 -0.928
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.2978828 -0.6328086 -0.1662572 0.5546343 4.0940868
##
## Number of Observations: 445
## Number of Groups: 17
##m4
m4 <- update(m3, correlation=corAR1(form = ~ 1 | ID))
summary(m4)$coef
## $fixed
## A.(Intercept) A.Read B.(Intercept) B.Read
## -2.74734459 7.22293656 0.01193214 0.04959430
##
## $random
## $random$ID
## A.(Intercept) B.(Intercept)
## S113 -4.085972e-09 -0.0270401937
## S110 6.172083e-07 -0.0337839074
## S101 -5.643789e-07 -0.0648744083
## S102 2.413142e-07 -0.0009086372
## S103 -2.609555e-08 -0.0540500740
## S104 -3.181928e-07 0.0525316734
## S108 -2.533298e-08 0.0778790059
## S116 -3.087107e-07 0.0486219917
## S106 8.287331e-08 -0.0727529894
## S105 1.920052e-07 0.0035211630
## S117 2.764930e-07 0.0494357469
## S115 -1.579409e-08 0.0296559402
## S107 -3.029171e-07 0.0029401752
## S114 3.284553e-08 0.0487637620
## S112 1.309247e-07 -0.0476460667
## S111 5.234170e-08 -0.0267635683
## S109 -6.049783e-08 0.0144703865
dta$yhat <- fitted(m3)
ggplot(dta, aes(Game, Moves), group = Read, coloer = Read ) +
geom_point(size = rel(.5), pch = 1, col = "skyblue") +
geom_line(aes(y = yhat, group = ID), col = "indianred") +
labs(x = "Number of game played (bouts)",
y = "Number of moves made before defeat (steps)") +
facet_wrap(~ ID, ncol = 6)+
theme(legend.position = "NONE")
這張圖顯示:除了前述之S101和S106與其他的ID不同外,S103的關聯程度亦不高。
#The End