[R-sig-ME] Request help translating LISREL code to R?

[Paul Johnson]

Hello again. I'm back, begging for help again.

I'm auditing a workshop on hierarchical models for psychologists and
they are using LISREL to estimate these models.  I volunteered to
translate the examples into R (either lme4 or nlme, as the case
requires).  Many of the examples are easily translate (some easy ones
pasted in at the end).

The ones that are trouble for me have more complicated covariance
structures among the random effects.

In LISREL, one can simply write down the covariance matrix to be
estimated, so, for example, if we write

1
2 3
0 0 6

it would mean we are setting the first 2 covariances in the last row
at 0 and all the rest are free.   Similarly, we can restrict estimates
to be the same by putting the same number in, like

1
0 3
0 0 3

Would force the last two diagonals to be equal.  Their variance
matrices look somewhat "uneven" because the matrix automatically
includes all random effects at a level.

This is data about children attitudes toward their parents, measured
at several time points.  The more I look at this, the less I
understand about the specific example, but I'm just replicating the
LISREL code, which is:

## Simchild_Lag1AutoCor.pr2
## OPTIONS OLS=YES CONVERGE=0.001000 MAXITER=10 OUTPUT=STANDARD;
 ## TITLE=Simulated Child Data;
 ## SY='simchild.psf';
 ## ID2=ID;
 ## RESPONSE=CLOSEDAD;
 ## FIXED=intcept GRADE_C;
 ## RANDOM2=intcept GRADE_C CONS1 CONS3 CONS4 CONS5 CONS6;
 ## RANDOM1=intcept;
 ## COV1PAT=0;
 ## COV2PAT=1
 ##         2 3
 ##         0 0 6
 ##         0 0 9 6
 ##         0 0 0 9 6
 ##         0 0 0 0 9 6
 ##         0 0 0 0 0 9 6;
 ## MISSING_DAT=-999999;

Here is the LISREL output:


 ITERATION NUMBER     6
                       +-----------------------+
                       | FIXED PART OF MODEL   |
                       +-----------------------+
------------------------------------------------------------------------------
  COEFFICIENTS             BETA-HAT      STD.ERR.      Z-VALUE       PR > |Z|
------------------------------------------------------------------------------
  intcept                 36.26858       0.11947     303.57744       0.00000
  GRADE_C                 -0.52397       0.03024     -17.32693       0.00000
                       +-----------------------+
                       |   -2 LOG-LIKELIHOOD   |
                       +-----------------------+
  DEVIANCE= -2*LOG(LIKELIHOOD) =  19383.33626094234
  NUMBER OF FREE PARAMETERS =     7
                       +-----------------------+
                       | RANDOM PART OF MODEL  |
                       +-----------------------+
------------------------------------------------------------------------------
  LEVEL 2                        TAU-HAT      STD.ERR.     Z-VALUE   PR > |Z|
------------------------------------------------------------------------------
  intcept /intcept               6.25166      0.70948      8.81167   0.00000
  GRADE_C /intcept               0.54434      0.14153      3.84619   0.00012
  GRADE_C /GRADE_C               0.12388      0.04831      2.56426   0.01034
  CONS1   /CONS1                 6.64178      0.31733     20.93050   0.00000
  CONS3   /CONS1                 0.25937      0.26284      0.98681   0.32374

        LEVEL 2 COVARIANCE MATRIX

               intcept    GRADE_C      CONS1      CONS3      CONS4      CONS5

 intcept       6.25166
 GRADE_C       0.54434    0.12388
 CONS1         0.00000    0.00000    6.64178
 CONS3         0.00000    0.00000    0.25937    6.64178
 CONS4         0.00000    0.00000    0.00000    0.25937    6.64178
 CONS5         0.00000    0.00000    0.00000    0.00000    0.25937    6.64178
 CONS6         0.00000    0.00000    0.00000    0.00000    0.00000    0.25937

                 CONS6
 CONS6         6.64178

        LEVEL 2 CORRELATION MATRIX

              intcept   GRADE_C     CONS1     CONS3     CONS4     CONS5
 intcept       1.0000
 GRADE_C       0.6186    1.0000
 CONS1         0.0000    0.0000    1.0000
 CONS3         0.0000    0.0000    0.0391    1.0000
 CONS4         0.0000    0.0000    0.0000    0.0391    1.0000
 CONS5         0.0000    0.0000    0.0000    0.0000    0.0391    1.0000
 CONS6         0.0000    0.0000    0.0000    0.0000    0.0000    0.0391


I'm trying to estimate that with lme.  (In their examples where lmer
can be used, this is easier for me).

In case you want the data, I uploaded it in spss format:

http://pj.freefaculty.org/R/simchild.sav

Here's how I'm getting started:

library(foreign)
dat <- read.spss("simchild.sav", to.data.frame=T)
oldnames <- colnames(dat)
colnames(dat) <-  tolower(oldnames)

## does not work-estimates not remotely same as LISREL
## library(nlme)
## mod7lme <- lme(closedad ~  grade_c, data=dat, random= ~ 1 + grade_c
| id,  correlation = corAR1(0.8, form= ~ 1 | id ))
## summary(mod7lme)

Thanks in advance, sorry if this is annoying.

-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas


Now, for some of the examples that DO work--lmer does match LISREL
output exactly



## Hox_Pop_Nested1.pr2
##  OPTIONS OLS=YES CONVERGE=0.001000 MAXITER=10 OUTPUT=STANDARD;
##  TITLE=Model 1 Nested Model Comparison;
##  SY='popular.psf';
##  ID2=SCHOOL;
##  RESPONSE=POPULAR;
##  FIXED=intcept SEX TEXP;
##  !boy=0 girl=1
##  RANDOM1=intcept;
##  RANDOM2=intcept;

mod2 <- lmer(popular ~ sex + texp +  (1| school), data=dat, REML=F)
summary(mod2)
anova(mod2)
## Matches LISREL Output



## SB_Add_Mean.PR2
##  OPTIONS OLS=YES CONVERGE=0.0000001000 MAXITER=20 OUTPUT=STANDARD;
##  TITLE=Uncentered Predictor Verbal IQ;
##  SY='SB.psf';
##  ID2=schoolnr;
##  RESPONSE=langpost;
##  FIXED=intcept IQV_CWC IQV_CM_C;
##  RANDOM1=intcept;
##  RANDOM2=intcept IQV_CWC;




mod4 <- lmer(langpost ~ IQV_CWC + IQV_CM_C + (1+ IQV_CWC | schoolnr),
data=dat, REML=F)
summary(mod4)
anova(mod4)
##Matches LISREL Output