fox & geese game
Szu-Yu Chen
04 Dec 2020
1 Data Management
load data
dtaHW2 <- read.table("C:/Users/ASUS/Desktop/data/fox_geese.txt", h=T)
head(dtaHW2)## ID Game Moves Read
## 1 S101 1 4 1.4
## 2 S101 2 7 1.4
## 3 S101 3 8 1.4
## 4 S101 4 3 1.4
## 5 S101 5 3 1.4
## 6 S101 6 3 1.42 Baisc graph
進攻次數與局數的關係(個別玩家),散佈圖
library(lattice)
xyplot(Moves ~ Game | ID , data=dtaHW2,
type=c("g", "p"), cex=.5,
xlab="Number of game played",
ylab="Number of moves made before defeat")
進攻次數與局數(所有玩家)，折線圖
library(ggplot2)
# individual profiles on one panel
ggplot(dtaHW2, aes(Game, Moves, group = ID)) +
geom_point(pch = 20)+
geom_line(alpha = .3) +
labs(x = "Number of game played", y = "Number of moves made before defeat")
進攻次數與局數的關係(個別玩家),折線圖
# coerce data into groupedData object
library(nlme)
dtaHW2_g <- groupedData(Moves ~ Game | ID, outer = ~ Read,
data = as.data.frame(dtaHW2),
labels = list(x = "Number of game played",
y = "Number of moves made before defeat"),
units = list(x = "(bouts)", y = "(steps)"))
# lattice plot by subject
plot(dtaHW2_g, pch = 1, aspect = .9)
3 non-linear Model
Define function as Yij = 1 + 19 / (1 + b0i exp(-b1iTimeij)) + εij, subject i (玩家) jth game played(局數) Yij (進攻次數) b0i (150, 15, 1.5) (A) b1i (0.1, 0.3, 0.5) (B)
定義函數mfx為上式右半邊(定義x, A, B 間的關係)
mfx <- function(x, A, B) { (1 + (19 /(1 + (A*exp(-B*x)))))}3.1 m0
m0以mfx來預測Moves(進攻次數)，定義A以150, B以0.1 為起始值
m0 <- nlsList(Moves ~ mfx(Game, A, B),
start = c(A = 150, B = 0.1),
data = dtaHW2_g)## Warning: 1 error caught in nls(model, data = data, control = controlvals, start
## = start): number of iterations exceeded maximum of 50Generate coefficients
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.1555140173.2 m1
Nonlinear mixed-effects model fit by maximum likelihood define A, B as uncorrelated parameters
m1 <- nlme(m0, random=pdDiag(A + B ~ 1))
m1## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Moves ~ mfx(Game, A, B)
## Data: dtaHW2_g
## Log-likelihood: -1245.006
## Fixed: list(A ~ 1, B ~ 1)
## A B
## 10.3571220 0.1093886
##
## Random effects:
## Formula: list(A ~ 1, B ~ 1)
## Level: ID
## Structure: Diagonal
## A B Residual
## StdDev: 2.628445 0.04992844 3.745158
##
## Number of Observations: 445
## Number of Groups: 17prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)
prm0## A B
## 10.357 0.1093.3 m2
in m2, update m1 with:
a constant variance function structure (grouped by ID) where A, B change to 8.57 and 0.096
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: dtaHW2_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: 173.4 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: dtaHW2_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: 173.5 m4
in m4, update m3 with:
correlation matrix as AR1
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.01447038653.6 m5
in m5, update m4 with:
only the matrix of A has diagnoal =1
m5 <- update(m4, random=pdDiag(A ~ 1))3.7 model comparison
Result showed m3 has better fit
anova(m1, m2, m3, m4, m5)## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 5 2500.012 2520.502 -1245.006
## m2 2 21 2456.174 2542.234 -1207.087 1 vs 2 75.83759 <.0001
## m3 3 23 2452.096 2546.352 -1203.048 2 vs 3 8.07780 0.0176
## m4 4 24 2423.659 2522.013 -1187.830 3 vs 4 30.43706 <.0001
## m5 5 23 2494.401 2588.656 -1224.200 4 vs 5 72.74149 <.00013.8 plots for each player(m3)
dtaHW2$yhat <- fitted(m3)
ggplot(dtaHW2, aes(Game, Moves), group = Read, color = 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")