r-crossover
Differences
This shows you the differences between two versions of the page.
| — |
r-crossover [2017/09/20 19:05] (current) |
||
|---|---|---|---|
| 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)[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)))) | ||
| + | </code> | ||
r-crossover.txt ยท Last modified: 2017/09/20 19:05 (external edit)
Page Tools
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 4.0 International
