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@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 ______________________________________________ R-help@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.
