I am trying to fit a simple linear mixed model with a intercept and
slope (fixed and random effects) where both are stratified by another
variable.

For example, suppose Xi is a subject specific binary variable and tj is
time for the jth observation, then the model I wish to fit is 

 

Yij=B0*Xi+B1tj*Xi+B2(1-Xi)+B2tj*(1-Xi)+bi0*Xi+bi1tjXi+bi2*(1-Xi)+bi3*tj*
(1-Xi), 

 

Where the random effects are block diagonal (i.,e. bi0 and bi1  are
correlated and bi2 and bi3 are correlated, but all other pairwise
combinations are 

Independent from each other). 

 

I wish to assign the correlation estimates to a function. When I do this
in SPLUS the third column produced in VarCorr provides the estimated
correlations exactly as I want them.

When I run the same code in R (which is what I want to use for my
simulations), the correlations are not provided. 

 

Here is the code that I run on SPLUS and R:

 

set.seed(119)

 

sigma0=0.75

sigmab0=1

sigmab1=1

 

 

ts=seq(1:5)

ID=rep(seq(1:3000),each=5)

 

 

for (ii in 1:1) {

print(ii)

 

X=c(rep(0,1500),rep(1,1500))

 

 

bi0=rnorm(3000)*sigmab0*(X+1)

bi10=rep(bi0,each=5)

 

 

bi1=rnorm(3000)*sigmab1*(X+1)

 

bi11=rep(bi1,each=5)

 

e1=rnorm(15000)*sigma0

 

tt=rep(seq(1:5),3000)

XX=rep(X,each=5)

 

y=1+XX+1*tt*(XX)+bi10+bi11*tt+e1

 

dum0=(XX==0)*tt

dum1=(XX==1)*tt

dum0b=(XX==0)*1

dum1b=(XX==1)*1

dat=cbind(ID,y,dum0,dum1,dum0b,dum1b,tt)

dat=data.frame(dat)

test=groupedData(y~-1+dum0b+dum1b+dum0+dum1|ID,data=dat)

ans=lme(y~-1+dum0b+dum1b+dum0+dum1,

random=pdBlocked(list(pdSymm(~dum0b+dum0-1),pdSymm(~dum1b+dum1-1))),

data=test)  

    

print(summary(ans))

print(VarCorr(ans))

 

 

}

 

The output from SPLUS is exactly what I want: 

 

Output of summary(ans):

Linear mixed-effects model fit by REML

 Data: test 

       AIC      BIC    logLik 

  54429.69 54513.46 -27203.84

 

Random effects:

 Composite Structure: Blocked

 

  Block 1: dum0b, dum0

 Formula:  ~ dum0b + dum0 - 1 | ID

 Structure: General positive-definite

         StdDev   Corr 

dum0b 0.9823633 dum0b 

 dum0 1.0086996 -0.012

 

  Block 2: dum1b, dum1

 Formula:  ~ dum1b + dum1 - 1 | ID

 Structure: General positive-definite

            StdDev  Corr 

   dum1b 1.9403494 dum1b

    dum1 2.0571202 0.006

Residual 0.7538521      

 

Fixed effects: y ~ -1 + dum0b + dum1b + dum0 + dum1 

          Value  Std.Error    DF   t-value p-value 

dum0b  1.023360 0.03255928  2998  31.43067  <.0001

dum1b  2.017929 0.05409916  2998  37.30055  <.0001

 dum0 -0.002711 0.02676196 11999  -0.10131  0.9193

 dum1  0.986042 0.05347007 11999  18.44100  <.0001

 

Standardized Within-Group Residuals:

      Min         Q1         Med        Q3      Max 

 -2.97078 -0.5454272 0.002605002 0.5435903 3.417185

 

 

Output of VarCorr:

 

Number of Observations: 15000

Number of Groups: 3000 

ID = pdBlocked(pdSymm(~ dum0b + dum0 - 1), pdSymm(~ dum1b + dum1 - 1)) 

          Variance    StdDev  Corr. 

   dum0b 0.9650377 0.9823633 dum0b 

    dum0 1.0174748 1.0086996 -0.012

   dum1b 3.7649557 1.9403494 dum1b 

    dum1 4.2317433 2.0571202 0.006 

Residual 0.5682930 0.7538521       

 

 

When I run this in R,  the output for VarCorr does not provide the Corr:

 

 

Number of Observations: 15000

Number of Groups: 3000 

ID = pdSymm(dum0b + dum0 - 1), pdSymm(dum1b + dum1 - 1) 

         Variance  StdDev   

dum0b    0.9785806 0.9892324

dum0     0.9915580 0.9957700

dum1b    3.9632302 1.9907863

dum1     3.9872415 1.9968078

Residual 0.5616301 0.7494198

 

The differences in the actual estimates are due to the fact that data
are generated within the program and the common seed generated a
different set of random variables 

For R and SPLUS.  

 

Would anyone know how to assign the correlation estimates to a function
in R?

 

 

Thanks!

 

Paul Albert

 

 

 

 

 

Paul S. Albert, Ph.D.

Biometric Research Branch

National Cancer Insitute

 


	[[alternative HTML version deleted]]


