Statistical aspects of two-way cross-over studies
[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
I ran into this presentation on Statistical aspects of two-way cross-over studies by Ing. Helmut Schütz (http://bebac.at). He presented some code and referred to the bear package. The bear package is menu driven, which is not my thing. I had to try and do that in R via other packages. The aim is to estimate if the treatments are bioequivalent and to estimate intra subject variation.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Data
This seems like a funny way to create data, but the presentation contains the data in these four lines. The code after that, putting it in a data.frame, is mine.
T1 <- c(28.39,33.36,24.29,28.61,59.49,42.30)
T2 <- c(49.42,36.78,34.81,45.54,28.23,25.71)
R1 <- c(39.86,32.75,34.97,45.44,27.87,24.26)
R2 <- c(35.44,33.40,24.65,31.77,65.29,37.01)
datain <- data.frame(
aval=c(T1,R2,R1,T2),
subjid=factor(c(
rep(c(1,4,6,7,9,12),2),
rep(c(2,3,5,8,10,11),2))),
seqa=rep(c(‘TR’,’RT’),each=12),
trta=rep(c(‘T’,’R’,’R’,’T’),each=6),
aperiod=factor(rep(c(1,2,1,2),each=6))
)
datain <- datain[order(datain$subjid,datain$aperiod),]
datain$logAval <- log(datain$aval)
Anova
With these balanced data Anova gives no problems
A1 <- aov(logAval ~ aperiod +seqa+ trta + Error(subjid),
data=datain)
summary(A1)
Error: subjid
Df Sum Sq Mean Sq F value Pr(>F)
seqa 1 0.0023 0.0023 0.014 0.907
Residuals 10 1.5943 0.1594
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
aperiod 1 0.02050 0.020501 3.784 0.0804 .
trta 1 0.00040 0.000397 0.073 0.7921
Residuals 10 0.05417 0.005417
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Extraction of parameters of interest.
The degrees of freedom are calculated outside the anova. This particular bit should be correct even with unbalanced data, but that is untested.
wi <- summary(A1)$`Error: Within`[[1]]
countn <- rowSums(table(datain$seqa,datain$subjid)>0)
sn <- sqrt(1/(2*countn[1]) + 1/(2*countn[2]))
SE <- sqrt(wi$`Mean Sq`[rownames(wi)=='Residuals'])*sn
td <- qt(.05,wi$`Df`[rownames(wi)=='Residuals'])
cat(‘Point estimate (%)’,100*exp(coef(A1)$Within[‘trtaT’]),’\n’,
‘90% interval (%) ‘ ,
100*exp(coef(A1)$Within[‘trtaT’]+td*SE),’to ‘,
100*exp(coef(A1)$Within[‘trtaT’]-td*SE),
‘\n’,
‘CVIntra (%) ‘,
100*sqrt(exp(wi$`Mean Sq`[rownames(wi)==’Residuals’])-1),
‘\n’)
Point estimate (%) 100.8168
90% interval (%) 95.47312 to 106.4596
CVIntra (%) 7.370138
nlme
nlme gave no problems. the parameters are the same.
library(nlme)
L1 <- lme(logAval ~ aperiod + seqa + trta,
random=~1|subjid,
data=datain)
summary(L1)
Linear mixed-effects model fit by REML
Data: datain
AIC BIC logLik
6.767892 12.74229 2.616054
Random effects:
Formula: ~1 | subjid
(Intercept) Residual
StdDev: 0.2775045 0.07360159
Fixed effects: logAval ~ aperiod + seqa + trta
Value Std.Error DF t-value p-value
(Intercept) 3.510257 0.11720776 10 29.949011 0.0000
aperiod2 0.058454 0.03004772 10 1.945360 0.0804
seqaTR 0.019577 0.16301059 10 0.120097 0.9068
trtaT 0.008135 0.03004772 10 0.270732 0.7921
Correlation:
(Intr) aperd2 seqaTR
aperiod2 -0.128
seqaTR -0.695 0.000
trtaT -0.128 0.000 0.000
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.21397837 -0.41862241 -0.05987476 0.37505780 1.30243785
Number of Observations: 24
Number of Groups: 12
Extraction of parameters of interest
One noticeable things here is that VarCorr() resulted in a character variable, so this is transformed by as.numerical.
coL1 <- fixef(L1)
SEL1 <- sqrt(vcov(L1)['trtaT','trtaT'])
dfL1 <- L1$fixDF$terms[names(L1$fixDF$terms)=='trta']
tdL1 <- qt(.05,dfL1)
cat(‘Point estimate (%)’,100*exp(coL1[‘trtaT’]),’\n’,
‘90% interval (%) ‘ ,
100*exp(coL1[‘trtaT’]+tdL1*SEL1),’to ‘,
100*exp(coL1[‘trtaT’]-tdL1*SEL1),
‘\n’,
‘CVIntra (%) ‘,
100*sqrt(exp(as.numeric(VarCorr(L1)[‘Residual’,’Variance’]))-1),
‘\n’)
Point estimate (%) 100.8168
90% interval (%) 95.47312 to 106.4596
CVIntra (%) 7.370138
MCMCglmm
Definitely not the standard or expected calculation, but the answers are really close to the expected values.M1 <- MCMCglmm(logAval ~ aperiod + seqa + trta,
random= ~ subjid ,
data=datain,family=’gaussian’,nitt=500000,thin=20,
burnin=50000,verbose=FALSE)
summary(M1)
Iterations = 50001:499981
Thinning interval = 20
Sample size = 22500
DIC: -41.6716
G-structure: ~subjid
post.mean l-95% CI u-95% CI eff.samp
subjid 0.09517 0.02514 0.1989 22500
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.006832 0.00187 0.01401 22500
Location effects: logAval ~ aperiod + seqa + trta
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 3.510768 3.245017 3.764446 22518 <4e-05 ***
aperiod2 0.058452 -0.008879 0.125118 22500 0.0838 .
seqaTR 0.018367 -0.326948 0.397998 22500 0.9148
trtaT 0.007888 -0.056976 0.077410 22500 0.7985
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Extraction of parameters of interest
Fairly straightforward, it is all in the samples. Notice the CV Intra is larger than expected, this was observed before in previous REML exercises. By and large, getting decent (point) estimates for variances is a major weakness for MCMCglmm.
M1est <- 100*exp(quantile(M1$Sol[,'trtaT'],c(.05,.5,.95)))
cat(‘Point estimate (%)’,M1est[2],’\n’,
‘90% interval (%) ‘ ,
M1est[1],’to ‘,M1est[3],
‘\n’,
‘CVIntra (%) ‘,
100*sqrt(exp(mean(M1$VCV[,’units’]))-1),
‘\n’)
Point estimate (%) 100.7712
90% interval (%) 95.45551 to 106.4557
CVIntra (%) 8.279958
lme4
lme4 was non-cooperating. I expected to use mcmcsamp(), but that function has disappeared in my version of lme4. Since I was a minor version behind on R I updated, which lead to updating Eclipse to from Juno to Kepler, which lead to updating Java while I was at it, but mcmcsamp() was gone. As a note, I had an old version of R/ & lme4 sitting around, the values out of mcmcsamp were quite different from expected. I am using version 0.999999-2 for these calculations.
library(lme4)
R1 <- lmer(logAval ~ aperiod + seqa + trta + (1|subjid),
data=datain)
summary(R1)
Linear mixed model fit by REML [‘lmerMod’]
Formula: logAval ~ aperiod + seqa + trta + (1 | subjid)
Data: datain
REML criterion at convergence: -5.2321
Random effects:
Groups Name Variance Std.Dev.
subjid (Intercept) 0.077009 0.2775
Residual 0.005417 0.0736
Number of obs: 24, groups: subjid, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.510257 0.117208 29.949
aperiod2 0.058454 0.030048 1.945
seqaTR 0.019577 0.163011 0.120
trtaT 0.008135 0.030048 0.271
Correlation of Fixed Effects:
(Intr) aperd2 seqaTR
aperiod2 -0.128
seqaTR -0.695 0.000
trtaT -0.128 0.000 0.000
Extraction of parameters of interest
The easy way out as detailed in http://glmm.wikidot.com/faq would be to impute the df value from another program, such as aov(). This would mean 10 df, and all the values would be equal to the presentation. I this blog I want to try something else. I wonder of this combination of simulate() and refit() is legit? The 90% interval is slightly narrower than the one in the presentation.
simeff <- function(m0) {
s <- simulate(m0)
fixef(refit(m0,s))[‘trtaT’]
}
ff <- replicate(1000,simeff(R1))
R1est <- 100*exp(quantile(ff,c(.05,.5,.95)))
cat(‘Point estimate (%)’,R1est[2],’\n’,
‘90% interval (%) ‘ ,
R1est[1],’to ‘,R1est[3],
‘\n’,
‘CVIntra (%) ‘,
100*sqrt(exp(sigma(R1)^2)-1),
‘\n’)
Point estimate (%) 100.7082
90% interval (%) 95.75779 to 105.7563
CVIntra (%) 7.370137
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.