[R] sem package and growth curves
Chuck Cleland
ccleland at optonline.net
Wed Mar 3 17:41:16 CET 2010
Dear John,
Thanks very much for your message. I should have looked at the help
page for sem() more closely. Thanks again for your excellent work on
the package.
Regards,
Chuck
On 3/3/2010 10:18 AM, John Fox wrote:
> 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.
>
> ______________________________________________
> 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
More information about the R-help
mailing list