Chuck and John, Thank you both for your help. I figured that my problem was trying to work through a new area for me, and trying to learn a new package for that area at the same time. You have both provided examples that clarified things that I either didn't understand about SEM in general or had overlooked in the documentation for the sem package.
Again, this has been very helpful, Dan Daniel Nordlund Bothell, WA USA > -----Original Message----- > From: John Fox [mailto:j...@mcmaster.ca] > Sent: Wednesday, March 03, 2010 7:19 AM > To: 'Chuck Cleland'; 'Daniel Nordlund' > Cc: 'r-help' > Subject: RE: [R] sem package and growth curves > > Dear Chuck and Daniel, > > First, thanks Chuck for fielding the question, which I didn't notice in > r-help. > > I can get solutions for models A, B, and C using the automatic start > values > along with the argument par.size="startvalues" to sem() (as recommended in > ?sem if there are convergence problems). For example, for Model A: > > -------- snip --------- > > > modA <- specify.model() > 1: I -> ALC1, NA, 1 > 2: I -> ALC2, NA, 1 > 3: I -> ALC3, NA, 1 > 4: S -> ALC1, NA, 0 > 5: S -> ALC2, NA, 0.75 > 6: S -> ALC3, NA, 1.75 > 7: UNIT -> I, Mi, NA > 8: UNIT -> S, Ms, NA > 9: I <-> I, Vi, NA > 10: S <-> S, Vs, NA > 11: I <-> S, Cis, NA > 12: ALC1 <-> ALC1, Vd1, NA > 13: ALC2 <-> ALC2, Vd2, NA > 14: ALC3 <-> ALC3, Vd3, NA > 15: > Read 14 records > > sem.modA <- sem(modA, alc2.modA.raw, 1122, fixed.x="UNIT", > par.size="startvalues", raw=TRUE) > > > > summary(sem.modA) > > Model fit to raw moment matrix. > > Model Chisquare = 0.048207 Df = 1 Pr(>Chisq) = 0.82621 > BIC = -6.9747 > > Normalized Residuals > Min. 1st Qu. Median Mean 3rd Qu. Max. > -0.04050 -0.03790 -0.01600 0.00603 0.03200 0.09620 > > Parameter Estimates > Estimate Std Error z value Pr(>|z|) > Mi 0.225625 0.0106901 21.1059 0.0000e+00 I <--- UNIT > Ms 0.035978 0.0073456 4.8979 9.6865e-07 S <--- UNIT > Vi 0.087039 0.0071035 12.2530 0.0000e+00 I <--> I > Vs 0.019764 0.0052178 3.7877 1.5205e-04 S <--> S > Cis -0.012476 0.0045780 -2.7251 6.4282e-03 S <--> I > Vd1 0.048428 0.0064146 7.5495 4.3743e-14 ALC1 <--> ALC1 > Vd2 0.075702 0.0044403 17.0488 0.0000e+00 ALC2 <--> ALC2 > Vd3 0.076698 0.0098901 7.7551 8.8818e-15 ALC3 <--> ALC3 > > Iterations = 57 > > -------- snip --------- > > Model D converges with the default setting of par.size: > > -------- snip --------- > > > alc2.modD.raw <- raw.moments(subset(alc2, > + select=c('PEER1','PEER2','PEER3','ALC1','ALC2','ALC3','UNIT'))) > > > > modD <- specify.model() > 1: Ia -> ALC1, NA, 1 > 2: Ia -> ALC2, NA, 1 > 3: Ia -> ALC3, NA, 1 > 4: Sa -> ALC1, NA, 0 > 5: Sa -> ALC2, NA, 0.75 > 6: Sa -> ALC3, NA, 1.75 > 7: UNIT -> Ia, Mia, NA > 8: UNIT -> Sa, Msa, NA > 9: Ip -> PEER1, NA, 1 > 10: Ip -> PEER2, NA, 1 > 11: Ip -> PEER3, NA, 1 > 12: Sp -> PEER1, NA, 0 > 13: Sp -> PEER2, NA, 0.75 > 14: Sp -> PEER3, NA, 1.75 > 15: Ip -> Ia, B1, NA > 16: Sp -> Ia, B2, NA > 17: Ip -> Sa, B3, NA > 18: Sp -> Sa, B4, NA > 19: UNIT -> Ip, Mip, NA > 20: UNIT -> Sp, Msp, NA > 21: Ia <-> Ia, Via, NA > 22: Sa <-> Sa, Vsa, NA > 23: Ia <-> Sa, Cisa, NA > 24: Ip <-> Ip, Vip, NA > 25: Sp <-> Sp, Vsp, NA > 26: Ip <-> Sp, Cisp, NA > 27: ALC1 <-> ALC1, Vd1, NA > 28: ALC2 <-> ALC2, Vd2, NA > 29: ALC3 <-> ALC3, Vd3, NA > 30: PEER1 <-> PEER1, Vd4, NA > 31: PEER2 <-> PEER2, Vd5, NA > 32: PEER3 <-> PEER3, Vd6, NA > 33: ALC1 <-> PEER1, Cd1, NA > 34: ALC2 <-> PEER2, Cd2, NA > 35: ALC3 <-> PEER3, Cd3, NA > 36: > Read 35 records > > sem.modD <- sem(modD, alc2.modD.raw, 1122, fixed.x=c("UNIT"), raw=TRUE) > > summary(sem.modD) > > Model fit to raw moment matrix. > > Model Chisquare = 11.557 Df = 4 Pr(>Chisq) = 0.020967 > BIC = -16.534 > > Normalized Residuals > Min. 1st Qu. Median Mean 3rd Qu. Max. > -0.91500 -0.39200 0.00105 0.09760 0.39900 1.61000 > > Parameter Estimates > Estimate Std Error z value Pr(>|z|) > Mia 0.0666214 0.0156727 4.25079 2.1302e-05 Ia <--- UNIT > Msa 0.0083040 0.0147616 0.56254 5.7375e-01 Sa <--- UNIT > B1 0.7985829 0.1028010 7.76824 7.9936e-15 Ia <--- Ip > B2 0.0804315 0.1840470 0.43702 6.6210e-01 Ia <--- Sp > B3 -0.1433386 0.0762547 -1.87973 6.0144e-02 Sa <--- Ip > B4 0.5766956 0.1938673 2.97469 2.9328e-03 Sa <--- Sp > Mip 0.1881743 0.0119530 15.74285 0.0000e+00 Ip <--- UNIT > Msp 0.0961698 0.0096929 9.92167 0.0000e+00 Sp <--- UNIT > Via 0.0421656 0.0074640 5.64920 1.6120e-08 Ia <--> Ia > Vsa 0.0092181 0.0054564 1.68941 9.1140e-02 Sa <--> Sa > Cisa -0.0063651 0.0051128 -1.24492 2.1316e-01 Sa <--> Ia > Vip 0.0696837 0.0103795 6.71357 1.8991e-11 Ip <--> Ip > Vsp 0.0284726 0.0089274 3.18936 1.4259e-03 Sp <--> Sp > Cisp 0.0011771 0.0071251 0.16521 8.6878e-01 Sp <--> Ip > Vd1 0.0480379 0.0063780 7.53177 4.9960e-14 ALC1 <--> ALC1 > Vd2 0.0762156 0.0044523 17.11821 0.0000e+00 ALC2 <--> ALC2 > Vd3 0.0762794 0.0097763 7.80249 5.9952e-15 ALC3 <--> ALC3 > Vd4 0.1057875 0.0108526 9.74770 0.0000e+00 PEER1 <--> PEER1 > Vd5 0.1712811 0.0087037 19.67904 0.0000e+00 PEER2 <--> PEER2 > Vd6 0.1289592 0.0177027 7.28471 3.2241e-13 PEER3 <--> PEER3 > Cd1 0.0109322 0.0061562 1.77578 7.5769e-02 PEER1 <--> ALC1 > Cd2 0.0339991 0.0046391 7.32874 2.3226e-13 PEER2 <--> ALC2 > Cd3 0.0374125 0.0101878 3.67229 2.4038e-04 PEER3 <--> ALC3 > > Iterations = 139 > > -------- snip --------- > > Regards, > John > > -------------------------------- > John Fox > Senator William McMaster > Professor of Social Statistics > Department of Sociology > McMaster University > Hamilton, Ontario, Canada > web: socserv.mcmaster.ca/jfox > > > > -----Original Message----- > > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] > On > > Behalf Of Chuck Cleland > > Sent: March-03-10 9:03 AM > > To: Daniel Nordlund > > Cc: 'r-help' > > Subject: Re: [R] sem package and growth curves > > > > 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. > ______________________________________________ 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.
