# Simulate a crossover design with the formula: # Response ~ 1 + Treatment + Order + Treatment:Order + (1 | Patient) # Fit simulated data with linear mixed effects model. Make # significance decision about treatment effect on the basis # of 95% confidence interval (i.e., significant if 95% CI fails # to include zero). # n - number of patients in each order group # sdW - within patient standard deviation # sdB - between patient standard deviation # beta - coefficient vector c(Intercept, Treatment, Order, Treatment:Order) simulate <- function(n, sdW=4, sdB=1, beta=c(8, 4, 0, 0), alpha=0.05) { require("lme4") Patient <- as.factor(rep(1:(2*n), rep(2, 2*n))) Treatment <- c(rep(c("Treatment1", "Treatment2"), n), rep(c("Treatment2", "Treatment1"), n)) Order <- rep(c("First", "Second"), 2*n) Data <- data.frame(Patient, Treatment, Order) CMat <- model.matrix(~ Treatment * Order + Patient, data=Data) Response <- CMat %*% c(beta, rnorm(2*n-1, 0, sdB)) + rnorm(4*n, 0, sdW) Data$Response <- Response Fit <- lmer(Response ~ (1 | Patient) + Treatment * Order, data=Data) Est <- fixef(Fit)[2] Ste <- sqrt(vcov(Fit)[2,2]) prod(Est + c(-1,1) * qnorm(1-alpha/2) * Ste) > 0 } # type I error for n=20 (result: 0.059) #mean(replicate(1000, simulate(n=20, beta=c(8, 0, 0, 0)))) # type I error for n=50 (result: 0.057) #mean(replicate(1000, simulate(n=50, beta=c(8, 0, 0, 0)))) # type I error for n=20 and order effect 2 (result: 0.062) #mean(replicate(1000, simulate(n=20, beta=c(8, 0, 2, 0)))) # type I error for n=50 and order effect 2 (result: 0.05) #mean(replicate(1000, simulate(n=50, beta=c(8, 0, 2, 0)))) # power for n=20 and treatment effect 4 (result: 0.869) #mean(replicate(1000, simulate(n=20, beta=c(8, 4, 0, 0)))) # power for n=50 and treatment effect 4 (result: 0.997) #mean(replicate(1000, simulate(n=50, beta=c(8, 4, 0, 0))))