[R] sem package and growth curves

[Chuck Cleland]

On 3/2/2010 1:43 AM, Daniel Nordlund wrote:
> I have been working through the book "Applied longitudinal data analysis: modeling change and event occurrence" by Judith D. Singer and John B. Willett.  I have been working examples using SAS and also using it as an opportunity for learning to use R for statistical analysis.
> 
> I ran into some difficulties in chapter 8 which deals with using structural equation modeling.  I have tried to use the sem package to replicate the problem solutions in chapter 8.  I am more familiar with RAM specifications than I am with structural equations (though I am a novice at both).  The solutions I have tried seem to be very sensitive to starting values (especially with more complex models).  I don't know if this is just my lack of knowledge in this area, or something else.
> 
> Has anyone worked out solutions to the Singer and Willett examples for Chapter 8 that they would be willing to share?  I would also be interested in other simple examples using sem and RAM specifications.  If anyone is interested, I would also be willing to share the R code I have written for other chapters in the Singer and Willett book.

Hi Dan,

  See below for my code for Models A-D in Chapter 8.  As you point out,
I find that this only works when good starting values are given.  I took
the starting values from the results given for another program (Mplus)
at the UCLA site for this text:

http://www.ats.ucla.edu/stat/examples/alda.htm

  I greatly appreciate John Fox's hard work on the sem package, but
since good starting values will generally not be available to applied
users I think the package is not as useful for these types of models as
it could be.  If anyone has approaches to specifying the models that are
less sensitive to starting values, or ways for less sophisticated users
to generate good starting values, please share.

Chuck

# Begin Code for Models A-D, Chapter 8, Singer & Willett (2003)

alc2 <-
read.table("http://www.ats.ucla.edu/stat/mplus/examples/alda/alcohol2.txt",
sep="\t", header=FALSE)

names(alc2) <- c('ID','FEMALE','ALC1','ALC2','ALC3','PEER1','PEER2','PEER3')

alc2$UNIT <- 1

library(sem)

alc2.modA.raw <- raw.moments(subset(alc2,
select=c('ALC1','ALC2','ALC3','UNIT')))

modA <- specify.model()
I -> ALC1, NA, 1
I -> ALC2, NA, 1
I -> ALC3, NA, 1
S -> ALC1, NA, 0
S -> ALC2, NA, 0.75
S -> ALC3, NA, 1.75
UNIT -> I, Mi, 0.226
UNIT -> S, Ms, 0.036
I <-> I, Vi, NA
S <-> S, Vs, NA
I <-> S, Cis, NA
ALC1 <-> ALC1, Vd1, 0.048
ALC2 <-> ALC2, Vd2, 0.076
ALC3 <-> ALC3, Vd3, 0.077

sem.modA <- sem(modA, alc2.modA.raw, 1122, fixed.x="UNIT", raw=TRUE)

summary(sem.modA)

alc2.modB.raw <- raw.moments(subset(alc2,
select=c('FEMALE','ALC1','ALC2','ALC3','UNIT')))

modB <- specify.model()
I -> ALC1, NA, 1
I -> ALC2, NA, 1
I -> ALC3, NA, 1
S -> ALC1, NA, 0
S -> ALC2, NA, 0.75
S -> ALC3, NA, 1.75
FEMALE -> I, B1, NA
FEMALE -> S, B2, NA
UNIT -> I, Mi, 0.226
UNIT -> S, Ms, 0.036
I <-> I, Vi, NA
S <-> S, Vs, NA
I <-> S, Cis, NA
ALC1 <-> ALC1, Vd1, 0.048
ALC2 <-> ALC2, Vd2, 0.076
ALC3 <-> ALC3, Vd3, 0.077

sem.modB <- sem(modB, alc2.modB.raw, 1122, fixed.x=c("FEMALE","UNIT"),
raw=TRUE)

summary(sem.modB)

alc2.modC.raw <- raw.moments(subset(alc2,
select=c('FEMALE','ALC1','ALC2','ALC3','UNIT')))

modC <- specify.model()
I -> ALC1, NA, 1
I -> ALC2, NA, 1
I -> ALC3, NA, 1
S -> ALC1, NA, 0
S -> ALC2, NA, 0.75
S -> ALC3, NA, 1.75
FEMALE -> I, B1, NA
FEMALE -> S, NA, 0
UNIT -> I, Mi, 0.226
UNIT -> S, Ms, 0.036
I <-> I, Vi, NA
S <-> S, Vs, NA
I <-> S, Cis, NA
ALC1 <-> ALC1, Vd1, 0.048
ALC2 <-> ALC2, Vd2, 0.076
ALC3 <-> ALC3, Vd3, 0.077

sem.modC <- sem(modC, alc2.modC.raw, 1122, fixed.x=c("FEMALE","UNIT"),
raw=TRUE)

summary(sem.modC)

alc2.modD.raw <- raw.moments(subset(alc2,
select=c('PEER1','PEER2','PEER3','ALC1','ALC2','ALC3','UNIT')))

modD <- specify.model()
Ia -> ALC1, NA, 1
Ia -> ALC2, NA, 1
Ia -> ALC3, NA, 1
Sa -> ALC1, NA, 0
Sa -> ALC2, NA, 0.75
Sa -> ALC3, NA, 1.75
UNIT -> Ia, Mia, 0.226
UNIT -> Sa, Msa, 0.036
Ip -> PEER1, NA, 1
Ip -> PEER2, NA, 1
Ip -> PEER3, NA, 1
Sp -> PEER1, NA, 0
Sp -> PEER2, NA, 0.75
Sp -> PEER3, NA, 1.75
Ip -> Ia, B1, 0.799
Sp -> Ia, B2, 0.080
Ip -> Sa, B3, -0.143
Sp -> Sa, B4, 0.577
UNIT -> Ip, Mip, 0.226
UNIT -> Sp, Msp, 0.036
Ia <-> Ia, Via, 0.042
Sa <-> Sa, Vsa, 0.009
Ia <-> Sa, Cisa, -0.006
Ip <-> Ip, Vip, 0.070
Sp <-> Sp, Vsp, 0.028
Ip <-> Sp, Cisp, 0.001
ALC1 <-> ALC1, Vd1, 0.048
ALC2 <-> ALC2, Vd2, 0.076
ALC3 <-> ALC3, Vd3, 0.077
PEER1 <-> PEER1, Vd4, 0.106
PEER2 <-> PEER2, Vd5, 0.171
PEER3 <-> PEER3, Vd6, 0.129
ALC1 <-> PEER1, Cd1, 0.011
ALC2 <-> PEER2, Cd2, 0.034
ALC3 <-> PEER3, Cd3, 0.037

sem.modD <- sem(modD, alc2.modD.raw, 1122, fixed.x=c("UNIT"), raw=TRUE)

summary(sem.modD)

> Thanks,
> 
> Dan
> 
> Daniel Nordlund
> Bothell, WA USA  
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894