r-crossover

# Differences

This shows you the differences between two versions of the page.

 — r-crossover [2017/09/20 19:05] (current) 2013/02/22 08:05 matt created 2013/02/22 08:05 matt created Line 1: Line 1: + <​code>​ + # 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) + 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)))) +
r-crossover.txt · Last modified: 2017/09/20 19:05 (external edit)

### Page Tools 