[R] sem package and growth curves
Daniel Nordlund
djnordlund at verizon.net
Wed Mar 3 18:22:31 CET 2010
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:jfox at 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-bounces at r-project.org [mailto:r-help-bounces at 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 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
> >
> > ______________________________________________
> > 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.
>
More information about the R-help
mailing list