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 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