[R] sem package and growth curves
John Fox
jfox at mcmaster.ca
Wed Mar 3 16:18:51 CET 2010
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