Matt's Notes

# 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))))