W12 Homework 2
Ching-Fang Wu
Mon Dec 07 12:02:28 2020
1 Introduction
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
2 Data management
# 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)3 Data Visualization
#個別受試者移動次數與遊戲之間的關係(依照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")
這張圖顯示出個體間的差異。
4 Modeling
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)
4.1 m0
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 50coef(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.0144703865dta$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