[R-sig-ME] Question About Cluster RCT analysis
Reinhold Kliegl
reinhold.kliegl at gmail.com
Thu Jul 22 22:08:09 CEST 2010
This is a really surprising result, but then I do not know anything
about t.test.cluster(). Anyway, I would recommend to artifiically
balance the clusters (either by randomly deleting or even simply
re-assigning cases) or by removing some outlier clusters completely.
Do lmer and geese reproduce the treatment effect for such data?
Reinhold Kliegl
On Thu, Jul 22, 2010 at 9:15 PM, Kevin E. Thorpe
<kevin.thorpe at utoronto.ca> wrote:
> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list