sem package and growth curves

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

sem package and growth curves

Daniel Nordlund-2
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.

Thanks,

Dan

Daniel Nordlund
Bothell, WA USA
 

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: sem package and growth curves

Chuck Cleland
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  
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: sem package and growth curves

Fox, John
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: [hidden email] [mailto:[hidden email]]
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
> >
> > ______________________________________________
> > [hidden email] 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
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: sem package and growth curves

Chuck Cleland
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: [hidden email] [mailto:[hidden email]]
> 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
>>>
>>> ______________________________________________
>>> [hidden email] 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
>>
>> ______________________________________________
>> [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: sem package and growth curves

Daniel Nordlund-2
In reply to this post by Fox, John
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:[hidden email]]
> 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: [hidden email] [mailto:[hidden email]]
> 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
> > >
> > > ______________________________________________
> > > [hidden email] 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
> >
> > ______________________________________________
> > [hidden email] 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.
>

______________________________________________
[hidden email] 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.