Supplemental Materials
An Introduction to Modeling Longitudinal Data With Generalized Additive Models: Applications to Single-Case Designs
by K. J. Sullivan et al., 2014, Psychological Methods
http://dx.doi.org/10.1037/met0000020
Note: Anything following ## is annotation of the code.
Example 1 Data: Tab Delimited Text File
x y z
1 12 0
2 3 0
3 9 0
4 10 0
5 4 0
9 8 1
11 2 1
16 11 1
18 7 1
23 8 1
25 30 1
30 9 1
32 14 1
37 11 1
39 13 1
44 22 1
46 12 1
51 21 1
53 21 1
58 2 1
In this data set, “x” is the trend/time variable, “y” is the outcome variable, and “z” is a treatment/baseline code, with 0=baseline and 1=treatment. The interaction term is calculated later.
l library(mgcv)
title <- read.table("C:/spss/kristynn sullivan cp/Example One 74 7 1/7472 time.txt",header=TRUE)
x <- title$x
x <- x-1
y <- title$y
z <- title$z
#trial <- title$trial
interaction <- function(x,ti,z) {
int <- (x-8)*z
}
int <- interaction(x,8,z)
gam1 <- gam(y~x+z+int, family=quasipoisson)
gam2 <- gam(y~x+z+s(int, bs="cr"), family=quasipoisson)
gam3 <- gam(y~s(x, bs="cr")+z+int, family=quasipoisson)
gam4 <- gam(y~s(x, bs = "cr")+z+s(int, bs = "cr"), family = quasipoisson)
summary(gam1)
gam.check(gam1)
gam1$coefficients
summary(gam2)
gam.check(gam2)
gam2$coefficients
summary(gam3)
gam.check(gam3)
gam3$coefficients
summary(gam4)
gam.check(gam4)
gam4$coefficients
anova(gam1,gam2, test="Chisq")
anova(gam1,gam3, test="Chisq")
anova(gam1,gam4, test="Chisq")
anova(gam2,gam4, test="Chisq")
anova(gam3,gam4, test="Chisq")
###for negative binomial estimation, substitute these lines for the quaipoisson lines above
gam1 <- gam(y~x+z+int, k = 4, family=negbin(1,link="log"), optimizer="perf")
gam2 <- gam(y~x+z+s(int, bs="cr", k = 4), family=negbin(1,link="log"), optimizer="perf")
gam3 <- gam(y~s(x, bs="cr", k = 4)+z+int, family=negbin(1,link="log"), optimizer="perf")
gam4 <- gam(y~s(x, bs = "cr", k = 4)+z+s(int, bs = "cr", k = 4), family=negbin(1,link="log"), optimizer="perf")
glm1 <- glm(y~x+z+int) ##models as normal distribution with trend
glm2 <- glm(y~z) ##models as normal dist no trend
summary(glm1)
glm1$coefficients
summary(glm2)
glm2$coefficients
plot(c(0,60), c(0,40), xlab="Day",ylab="Number of Peer Requests",type="n")
points(x,y,cex=1.25)
abline(v=7.5)
points(x,gam1$fit,type="o",pch=20)
points(x,gam2$fit,type="b",pch=10,lty=2)
points(x,glm2$fit,type="b",pch=4,lty=3)
Example 2
x y z trial
1 .09 0 11
2 .09 0 11
3 .09 0 11
4 .09 0 11
5 .09 0 11
8 .09 0 11
9 .09 0 11
10 .09 0 11
11 .09 0 11
12 .18 0 11
15 .18 0 11
16 .18 0 11
17 .09 1 11
18 .09 1 11
19 .09 1 11
22 .27 1 11
23 .09 1 11
24 .18 1 11
25 .09 1 11
26 .18 1 11
29 .18 1 11
30 .18 1 11
31 .18 1 11
32 .36 1 11
33 .27 1 11
36 .27 1 11
37 .27 1 11
38 .45 1 11
39 .45 1 11
40 .54 1 11
43 .54 1 11
44 .54 1 11
45 .63 1 11
46 .63 1 11
47 .72 1 11
50 .72 1 11
51 .81 1 11
52 .90 1 11
54 .90 1 11
57 .90 1 11
58 .81 1 11
60 .90 1 11
64 .90 1 11
65 .90 1 11
68 .90 1 11
71 .90 1 11
library(mgcv)
title <- read.table("C:/spss/kristynn sullivan cp/Example Three 105 1 2 now two/10512 time.txt",header=TRUE)
x <- title$x
x <- x-1
y <- title$y
z <- title$z
trial <- title$trial
interaction <- function(x,ti,z) {
int <- (x-16)*z
}
int <- interaction(x,16,z)
gam1 <- gam(y~x+z+int, family=quasibinomial, weights=trial)
gam2 <- gam(y~x+z+s(int, bs="cr"), family=quasibinomial, weights=trial)
gam3 <- gam(y~s(x, bs="cr")+z+int, family=quasibinomial, weights=trial)
gam4 <- gam(y~s(x, bs = "cr")+z+s(int, bs = "cr"), family=quasibinomial, weights=trial)
summary(gam1)
gam.check(gam1)
summary(gam2)
summary(gam3)
summary(gam4)
anova(gam1,gam2, test="Chisq")
anova(gam1,gam3, test="Chisq")
anova(gam1,gam4, test="Chisq")
anova(gam2,gam4, test="Chisq")
anova(gam3,gam4, test="Chisq")
glm1 <- glm(y~x+z+int)
glm2 <- glm(y~z)
summary(glm1)
summary(glm2)
gamm2.105.1.2 <- gamm(y~x+z+s(int,bs="cr"),family=quasipoisson, correlation=corAR1(form=~1|x))
gamm2.105.1.2
plot(c(0,75), c(0,1), xlab="Day",ylab="Percentage of Social Play Skills",type="n")
points(x,y,cex=1.25)
abline(v=15.5)
points(x,gam1$fit,type="o",pch=20)
points(x,gam2$fit,type="b",pch=10,lty=2)
points(x,glm2$fit,type="b",pch=4,lty=3)
Example 3
x y z
1 32 0
2 16 0
3 17 0
13 22 0
21 14 0
25 12 0
28 12 0
29 16 0
30 0 1
31 3 1
32 16 1
33 11 1
36 1 1
37 0 1
39 7 1
40 10 1
library(mgcv)
title <- read.table("C:/spss/kristynn sullivan cp/Example Four 80 1 5 now three/8015 session.txt",header=TRUE)
x <- title$x
x <- x-1
y <- title$y
z <- title$z
#trial <- title$trial
interaction <- function(x,ti,z) {
int <- (x-29)*z
}
int <- interaction(x,29,z)
gam1 <- gam(y~x+z+int, k = 4, family=quasipoisson,)
gam2 <- gam(y~x+z+s(int, bs="cr", k = 4), family=quasipoisson)
gam3 <- gam(y~s(x, bs="cr", k = 4)+z+int, family=quasipoisson)
gam4 <- gam(y~s(x, bs = "cr", k = 4)+z+s(int, bs = "cr", k = 4), family = quasipoisson)
summary(gam1)
gam1$coefficients
summary(gam2)
gam2$coefficients
summary(gam3)
gam3$coefficients
summary(gam4)
gam4$coefficients
anova(gam1,gam2, test="Chisq")
anova(gam1,gam3, test="Chisq")
anova(gam1,gam4, test="Chisq")
anova(gam2,gam4, test="Chisq")
anova(gam3,gam4, test="Chisq")
###for negative binomial estimation, substitute these lines for the quaipoisson lines above
gam1 <- gam(y~x+z+int, k = 4, family=negbin(1,link="log"), optimizer="perf")
gam2 <- gam(y~x+z+s(int, bs="cr", k = 4), family=negbin(1,link="log"), optimizer="perf")
gam3 <- gam(y~s(x, bs="cr", k = 4)+z+int, family=negbin(1,link="log"), optimizer="perf")
gam4 <- gam(y~s(x, bs = "cr", k = 4)+z+s(int, bs = "cr", k = 4), family=negbin(1,link="log"), optimizer="perf")
glm1 <- glm(y~x+z+int) ##models as normal distribution with trend
glm2 <- glm(y~z) ##models as normal dist no trend
summary(glm1)
glm1$coefficients
summary(glm2)
glm2$coefficients
plot(c(0,46), c(0,35), xlab="Session",ylab="Responses",type="n")
points(x,y,cex=1.25)
abline(v=28.5)
points(x,gam1$fit,type="o",pch=20)
points(x,gam4$fit,type="b",pch=10,lty=2)
points(x,glm2$fit,type="b",pch=4,lty=3)