[R-sig-ME] Question About Cluster RCT analysis

[Kevin E. Thorpe]

Hello.

I'm in the process of analyzing a cluster RCT with a continuous outcome
and am comparing some methods to help aid my understanding.  I have
compared the results from t.test.cluster in Hmisc to results using lmer
in lme4 and geese in geepack.  My main question concerns the effect size
(follows the output from all three analyses).  Here is the output from 
t.test.cluster.

> with(fim, t.test.cluster(FIM_TotalScore,Hospital_Code,Group))
                                    Control    Intervention
N                                   882        921
Clusters                            10         10
Mean                                98.00680   99.79045
SS among clusters within groups     36569.40   47721.81
SS within clusters within groups    411844.6   445172.7
MS among clusters within groups     4682.845
d.f.                                18
MS within clusters within groups    480.6603
d.f.                                1783
Na                                  82.28852
Intracluster correlation            0.09603894
Variance Correction Factor          12.63857   12.17243
Variance of effect                  13.24026
Variance without cluster adjustment 1.066856
Design Effect                       12.41054
Effect (Difference in Means)        1.783642
S.E. of Effect                      3.638716
0.95 Confidence limits              -5.348111   8.915396
Z Statistic                         0.4901845
2-sided P Value                     0.6240033

Now, lmer.

> lmer(FIM_TotalScore~Group+(1|Hospital_Code),data=fim)
Linear mixed model fit by REML
Formula: FIM_TotalScore ~ Group + (1 | Hospital_Code)
   Data: fim
   AIC   BIC logLik deviance REMLdev
16292 16314  -8142    16291   16284
Random effects:
Groups        Name        Variance Std.Dev.
Hospital_Code (Intercept)  46.279   6.8029
Residual                  480.657  21.9239
Number of obs: 1803, groups: Hospital_Code, 20

Fixed effects:
                  Estimate Std. Error t value
(Intercept)        98.7553     2.3091   42.77
GroupIntervention   0.5398     3.2653    0.17

Correlation of Fixed Effects:
            (Intr)
GrpIntrvntn -0.707
> 46.279/(46.279+480.657)  # icc

Finally, geese.

> 
summary(geese(FIM_TotalScore~Group,id=Hospital_Code,data=fim,corstr="exchangeable"))

Call:
geese(formula = FIM_TotalScore ~ Group, id = Hospital_Code, data = fim,
    corstr = "exchangeable")

Mean Model:
Mean Link:                 identity
Variance to Mean Relation: gaussian

Coefficients:
                    estimate   san.se         wald         p
(Intercept)       98.7547493 2.052763 2.314399e+03 0.0000000
GroupIntervention  0.5472461 3.195752 2.932373e-02 0.8640337

Scale Model:
Scale Link:                identity

Estimated Scale Parameters:
            estimate   san.se     wald p
(Intercept) 522.4746 50.14271 108.5712 0

Correlation Model:
Correlation Structure:     exchangeable
Correlation Link:          identity

Estimated Correlation Parameters:
       estimate     san.se     wald           p
alpha 0.0828722 0.02078484 15.89734 6.68725e-05

Returned Error Value:    0
Number of clusters:   20   Maximum cluster size: 223

As you can see, the effect size in t.test.cluser is 1.783642 which is
the difference in the means, which makes sense to me.  I would have
expected the estimate of the fixed group effect in lmer and geese to be
similar to this, which they are not, although they are similar to each
other.  lmer = 0.5398 and geese = 0.5472461.  The intercepts in both
cases are very close to the control group mean.  The icc estimates are
close, lmer = 0.0878266 (based on the random effect variances),
geese = 0.0828722 (the correlation in the exchangeable structure),
t.test.cluster = 0.09603894.

My contrasts are not weird.

> options("contrasts")
$contrasts
        unordered           ordered
"contr.treatment"      "contr.poly"

> contrasts(fim$Group)
             Intervention
Control                 0
Intervention            1

I'm probably being dense today, but can you offer any explanation for
the difference even if that involves exposing my denseness?

 > sessionInfo()
R version 2.10.1 Patched (2009-12-29 r50852)
i686-pc-linux-gnu

locale:
  [1] LC_CTYPE=en_US       LC_NUMERIC=C         LC_TIME=en_US
  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=en_US
  [7] LC_PAPER=en_US       LC_NAME=C            LC_ADDRESS=C
[10] LC_TELEPHONE=C       LC_MEASUREMENT=en_US LC_IDENTIFICATION=C

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] geepack_1.0-17     doBy_4.0.5         lme4_0.999375-33 
Matrix_0.999375-38
[5] lattice_0.18-3     Hmisc_3.7-0        survival_2.35-8

loaded via a namespace (and not attached):
[1] cluster_1.12.3 grid_2.10.1    nlme_3.1-96    tools_2.10.1


Kevin

-- 
Kevin E. Thorpe
Biostatistician/Trialist, Knowledge Translation Program
Assistant Professor, Dalla Lana School of Public Health
University of Toronto
email: kevin.thorpe at utoronto.ca  Tel: 416.864.5776  Fax: 416.864.3016