Dear R users,


I am a medical researcher at the University of Geneva undertaking
neuroscience research. I am writing to ask for some help with some feedback
from reviewers following submission of a paper leveraging the lme4 package
in R. The objective of the study was to understand the effect of congenital
heart disease requiring intervention on long term neurodevelopment.


First, a description of the data set.

Combination of data gathered longitudinally in France, the UK and
Switzerland.

Patients were collected as a part of the longitudinal cohort recruited in
Geneva since 2001.
80 children with DiGeorge syndrome
Measuring MRI brain volumes
Inclusion of repeated measurements for 44 out of the 80 patients, with up
to 4 scans per subject and providing a total of 164 scans for the study
Additional data on gender, scanning equipment and age

Following optimisation based on total brain volume I developed the
following formulas for each trajectory to be applied to each cortical and
hippocampal region before being re-run:

lm.Intercept<-lmer(Data~Group_CHD+Gender+Type_Scan+(1|SUBJ_ID),hippocampus,REML=TRUE)
lm.Linear<-lmer(Data~Mean_Age*Group_CHD+Mean_Age*Gender+Mean_Age*Type_Scan+(1|SUBJ_ID),hippocampus,REML=TRUE)
lm.Quadratic<-lmer(Data~Mean_Age*Group_CHD+Mean_Age*Gender+Mean_Age*Type_Scan+I(Mean_Age^2)*Group_CHD+I(Mean_Age^2)*Gender+I(Mean_Age^2)*Type_Scan+(1|SUBJ_ID),hippocampus,REML=TRUE)
lm.Cubic<-lmer(Data~Mean_Age*Group_CHD+Mean_Age*Gender+Mean_Age*Type_Scan+I(Mean_Age^2)*Group_CHD+I(Mean_Age^2)*Gender+I(Mean_Age^2)*Type_Scan+I(Mean_Age^3)*Group_CHD+I(Mean_Age^3)*Gender+I(Mean_Age^3)*Type_Scan+(1|SUBJ_ID),hippocampus,REML=TRUE)

I developed an algorithm to select a growth trajectory based on minimum AIC
for intercept, linear, quadratic and cubic mixed models, applying it to 34
cortical regions of the brain and 10 regions of the hippocampus in each
hemisphere. I ran the model with REML.

where fixed effects:
Data - Response vector of the volume in question
Mean_Age - Continuous variable of age centred on mean for the dataset (to
avoid autocorrelation)
Group_CHD - Binary value of a history of congenital heart disease
Gender - Binary value of gender
Type_Scan - Binary value of 1.5T or 3T studies

and one random effect of SUBJ_ID to account of inter-subject differences in
the intercept.

I then used pvals.fnc to extract p values for the difference in intercept
and ran Monte Carlo simulations to extract a "combined" p value of all
coefficients combining to form the "slope". For instance, for the cubic
trajectory, I used linear, quadratic and cubic terms together in
calculating a p-value with the following formula:

mcmcpvalue <- function(samp)
{
   ## elementary version that creates an empirical p-value for the
   ## hypothesis that the columns of samp have mean zero versus a
   ## general multivariate distribution with elliptical contours.
   ## differences from the mean standardized by the observed
   ## variance-covariance factor
   std <- backsolve(chol(var(samp)),
                    cbind(0, t(samp)) - colMeans(samp),
                    transpose = TRUE)
   sqdist <- colSums(std * std)
   sum(sqdist[-1] > sqdist[1])/nrow(samp)
}

I did this in order to investigate the absolute change in brain volume as
well as the way the brain developed over time .i.e. Mean_Age:Group_CHD,
I(Mean_Age^2):Group_CHD and I(Mean_Age)^3:Group_CHD.

The results showed for 17 cortical regions a significant difference in the
intercept term but not for any involved in the slopes of the trajectories,

The paper was accepted with revisions which is where we are now. The
reviewers correctly questioned to what extent we had sufficient power to
come to the conclusions we did. I do not want to perform a post-hoc power
analysis because from what I have read and seen on R discussions it is
discouraged. I also read the following papers supporting this decision.

Hoenig, J. M. & Heisey, D. M. The abuse of power: the pervasive fallacy of
power calculations for data analysis Amer Statistician, 2001, 55, 19-24
M Levine and MHH Ensom, Post Hoc Power Analysis: An Idea Whose Time Has
Passed <http://www.ncbi.nlm.nih.gov/pubmed/11310512>, Pharmacotherapy
21:405-409, 2001
Baayen, R. H.; Davidson, D. J. & Bates, D. M. Mixed-effects modeling with
crossed random effects for subjects and items J Mem Lang, 2008, 59, 390-412

I did however read of a simulation based sample size analysis which which
power can be extracted. This is something I pursued so that I had something
to report in the way of power. I used the following R formula to calculate
"power", based somewhat on a model for lme4 originally proposed by Greg
Snow (found here -
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q1/001790.html)

pboot <- function(m0,m1) {
  s0 <- simulate(m1)
  r0 <- refit(m0,s0)
  r1 <- refit(m1,s0)
  anova(r0,r1)[2,7]
}

Result<-replicate(1000,pboot(lm.Cubic3,lm.Cubic))

hist(Result)
mean(Result<=0.05)
prop.test(sum(Result<=0.05), 1000)

The formula "worked" without error with a suitable histogram of results.

Secondly, the reviewers requested confidence intervals where there was
significance in shape, but no following correction (Bonferroni). I used
HPD95upper and HPD95lower results from pvals.fnc in order to construct
confidence curves manually and superimpose them on manually constructed
trajectories and the raw data.

I am writing on here to ask if you agree with my approach in both
calculating a power number and constructing confidence intervals. If not,
could you suggest an alternative that would work with my unbalanced mixed
model (random intercept for subject and fixed effects for gender, scanning
equipment and history of CHD, with scans per subject ranging from 1-4 with
different time intervals).

If any of you are happy with the approach and want to see the algorithm or
code, I am happy to send!


Best,

Daniel Fountain
Researcher, University of Geneva
Medical Student, University of Cambridge
dmf34@cam.ac.uk

	[[alternative HTML version deleted]]

