# [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:
> 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:
> 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
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 <-
>
>
> 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
> 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