[R] Random intercept model with time-dependent covariates, results different from SAS

Keith Wong keithw at med.usyd.edu.au
Sun Jul 4 11:21:32 CEST 2004


Thank you for the very prompt response. I only included a small part of the
output to make the message brief. I'm sorry it did not provide enough detail to
answer my question. I have appended the summary() and anova() outputs to the
two models I fitted in R.

Quoting Prof Brian Ripley <ripley at stats.ox.ac.uk>:

> Looking at the significance of a main effect (group) in the presence of an
> interaction (time:group) is hard to interpret, and in your case is I think
> not even interesting.  (The `main effect' probably represents difference
> in intercept for the time effect, that is the group difference at the last
> time.  But see the next para.)  Note that the two systems are returning
> different denominator dfs.


I take your point that the main effect is probably not interesting in the
presence of an interaction. I was checking the results for consistency to see
if I was doing the right thing. I was not 100% sure that the SAS code was in
itself correct. 
 
> At this point you have not told us enough.  My guess is that you have
> complete balance with the same number of subjects in each group.  In that
> case the `group' effect is in the between-subjects stratum (as defined for
> the use of Error in aov, which you could also do), and thus R's 11 df
> would be right (rather than 44, without W and Z).  Without balance Type
> III tests get much harder to interpret and the `group' effect would appear
> in two strata and there is no simple F test in the classical theory.  So
> further guessing, SAS may have failed to detect balance and so used the
> wrong test.

I had not appreciated the need for balance: in actual fact, one group has 5
subjects and the other 7. Will this be a problem? Would the R analysis still be
valid in that case?

 
> The time-dependent covariates muddy the issue more, and I looked mainly at 
> the analyses without them.  Again, a crucial fact is not here: do the 
> covariates depend on the subjects as well?

Yes the covariates are measures of blood pressure and pulse, and they depend on
the subjects as well.

> The good news is that the results _are_ similar.  You do have different
> time behaviour in the two groups.  So stop worrying about tests of
> uninteresting hypotheses and concentrate of summarizing that difference.
> 
> -- 
> Brian D. Ripley,                  ripley at stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595


Thank you. I was concerned that one or both methods were incorrect given the
results were inconsistent. Perhaps reassuringly, the parameter estimates for
the fixed effects in both SAS and R were the same.

Is the model specification OK for the model with just time, group and their
interaction?
Is the model specification with the 2 time dependent covariates appropriate?

Once again, I'm very grateful for the time you've taken to answer my questions.

Keith

[Output from the 2 models fitted in R follows]

> g1 = lme(Y ~ time + group + time:group, random = ~ 1 | id, data = datamod)

> anova(g1)
            numDF denDF   F-value p-value
(Intercept)     1    44  3.387117  0.0725
time            4    44 10.620547  <.0001
group           1    11  0.508092  0.4908
time:group      4    44  3.961726  0.0079
> summary(g1)
Linear mixed-effects model fit by REML
 Data: datamod 
       AIC      BIC    logLik
  372.4328 396.5208 -174.2164

Random effects:
 Formula: ~1 | id
        (Intercept) Residual
StdDev:    11.05975 3.228684

Fixed effects: Y ~ time + group + time:group 
              Value Std.Error DF   t-value p-value
(Intercept)   8.250  4.073428 44  2.025321  0.0489
time1        -0.250  1.614342 44 -0.154862  0.8776
time2        -8.125  1.614342 44 -5.033011  0.0000
time3        -8.875  1.614342 44 -5.497596  0.0000
time4        -4.250  1.614342 44 -2.632652  0.0116
group1        2.126  6.568205 11  0.323681  0.7523
time1:group1 -2.734  2.603048 44 -1.050307  0.2993
time2:group1  5.583  2.603048 44  2.144793  0.0375
time3:group1  5.549  2.603048 44  2.131732  0.0387
time4:group1  3.634  2.603048 44  1.396056  0.1697
 Correlation: 
             (Intr) time1  time2  time3  time4  group1 tm1:g1 tm2:g1 tm3:g1
time1        -0.198                                                        
time2        -0.198  0.500                                                 
time3        -0.198  0.500  0.500                                          
time4        -0.198  0.500  0.500  0.500                                   
group1       -0.620  0.123  0.123  0.123  0.123                            
time1:group1  0.123 -0.620 -0.310 -0.310 -0.310 -0.198                     
time2:group1  0.123 -0.310 -0.620 -0.310 -0.310 -0.198  0.500              
time3:group1  0.123 -0.310 -0.310 -0.620 -0.310 -0.198  0.500  0.500       
time4:group1  0.123 -0.310 -0.310 -0.310 -0.620 -0.198  0.500  0.500  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.63416413 -0.42033405  0.03577472  0.46164486  1.74068368 

Number of Observations: 65
Number of Groups: 13 


> g2 = lme(Y ~ time + group + time:group + W + Z, random = ~ 1 | id, data =
datamod)

> anova(g2)
            numDF denDF  F-value p-value
(Intercept)     1    42  5.54545  0.0233
time            4    42 16.41069  <.0001
group           1    11  0.83186  0.3813
W               1    42  0.07555  0.7848
Z               1    42 45.23577  <.0001
time:group      4    42  3.04313  0.0273

> summary(g2)
Linear mixed-effects model fit by REML
 Data: datamod 
       AIC      BIC    logLik
  355.2404 382.8245 -163.6202

Random effects:
 Formula: ~1 | id
        (Intercept) Residual
StdDev:    8.639157 2.597380

Fixed effects: Y ~ time + group + time:group + W + Z 
                 Value Std.Error DF   t-value p-value
(Intercept)  10.056433  9.583658 42  1.049331  0.3000
time1         0.209668  1.301306 42  0.161121  0.8728
time2         4.111435  2.556420 42  1.608278  0.1153
time3         0.423056  2.077066 42  0.203679  0.8396
time4        -3.976417  1.300572 42 -3.057437  0.0039
group1        4.677706  5.162006 11  0.906180  0.3843
W             0.377142  0.127146 42  2.966212  0.0050
Z            -0.531895  0.093276 42 -5.702395  0.0000
time1:group1 -0.845857  2.126289 42 -0.397809  0.6928
time2:group1 -5.145361  2.962470 42 -1.736848  0.0897
time3:group1 -3.261241  2.597008 42 -1.255769  0.2161
time4:group1  4.153245  2.096587 42  1.980956  0.0542
 Correlation: 
             (Intr) time1  time2  time3  time4  group1 W      Z      tm1:g1
tm2:g1
time1        -0.051                                                             
 
time2         0.199  0.308                                                      
 
time3         0.023  0.361  0.817                                               
 
time4        -0.029  0.501  0.293  0.342                                        
 
group1       -0.202  0.131  0.136  0.146  0.129                                 
 
W            -0.790  0.019  0.243  0.366 -0.015  0.044                          
 
Z            -0.146 -0.063 -0.853 -0.779 -0.041 -0.086 -0.409                   
 
time1:group1 -0.028 -0.601 -0.043 -0.074 -0.302 -0.187  0.147 -0.144            
 
time2:group1 -0.293 -0.262 -0.818 -0.642 -0.255 -0.198 -0.051  0.665  0.276     
 
time3:group1 -0.016 -0.286 -0.626 -0.774 -0.273 -0.214 -0.277  0.590  0.308 
0.668
time4:group1  0.065 -0.306 -0.116 -0.159 -0.616 -0.199  0.002 -0.046  0.497 
0.318
             tm3:g1
time1              
time2              
time3              
time4              
group1             
W                  
Z                  
time1:group1       
time2:group1       
time3:group1       
time4:group1  0.376

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.11181231 -0.43210237  0.04949838  0.32444580  2.77710590 

Number of Observations: 65
Number of Groups: 13 
>




More information about the R-help mailing list