\documentclass[12pt]{article} \usepackage[margin=0.5in]{geometry} \renewcommand{\familydefault}{\sfdefault} \begin{document} \SweaveOpts{concordance=FALSE} \begin{center} {\Large Collaborative Grant Writing Exercise}\\\vspace{0.1in} BIOS 352 \\\vspace{0.1in} \today \end{center} \noindent Consider a research proposal that aims to test whether the pharmacokinetics of a new drug are systematically different in women versus men. In this study, 10 women and 10 men will be recruited in the PI's clinic. Every patient that visits the clinic will be approached for recruitment until 20 total participants have given consent. During their visit, each participant will recieve a single fixed dose of the new drug via intravenous bolus, and will be monitored for 20 minutes afterward. Drug concentration in a body tissue will be measured non-invasively at baseline (0 minutes) and every 2 minutes, up to 20 minutes following drug administration. The primary hypothesis is that the new drug will be eliminated faster in women than in men.\\ \noindent Drug concentration will be modeled using a solution to the two-compartment pharmacokinetic model, as discussed in class. All three rate parameters (i.e., $k_{12}$, $k_{21}$, and $k_{10}$) will be estimated for both males and females. The difference in the drug elimination rate estimates ($k_{10}$) will be used to test the primary hypothesis. Develop and write a statistical analysis plan and power analysis that addresses the primary hypothesis. For the purpose of this assignment, it is not necessary to implement either the statistical analysis plan or power analysis. Assume that the power analysis has already been completed, and that there is 80\% power to detect a 0.1 unit difference in the drug elimination rate. However, be sure to describe how the power analysis was (would have been) conducted. There should be sufficient detail that a knowledgeable reader could reproduce your work. The write-up should be no more than one page in length. Use this document as a page template (i.e., sans-serif, 12pt font, except mathematical notation; 0.5in margins).\\ \hrule\vspace{0.3cm} \noindent The remainder of this document gives an implementation of the "mixed-effects" approach described in class. It is not needed to complete the assignment, but may be helpful in thinking about the statistical approach, and the details that may be necessary for an outsider to implement the statistical plan. The following code chuck sets up functions to simulate data from the nonlinear bi-exponential mixed effects model, and fit the model using the {\tt nlme} function. \scriptsize <<>>= set.seed(33) # bi-exponential tissue concentration Ctissue <- function(t, k12, k21, k10, Vb=1, Vt=1, D=1000) { alp <- (k12+k21+k10-sqrt((k12+k21+k10)^2-4*k21*k10)) bet <- (k12+k21+k10+sqrt((k12+k21+k10)^2-4*k21*k10)) D/Vt*k12/(alp-bet)*(exp(-bet*t)-exp(-alp*t)) } # bi-exponential blood concentration CBlood <- function(t, k12, k21, k10, Vb=1, Vt=1, D=1000) { alp <- (k12+k21+k10-sqrt((k12+k21+k10)^2-4*k21*k10)) bet <- (k12+k21+k10+sqrt((k12+k21+k10)^2-4*k21*k10)) D/Vb*(k21-alp)/(bet-alp)*exp(-alp*t)+(k21-bet)/(alp-bet)*exp(-bet*t) } # plot example concentration curves (bi-exponential) plotExampleCbCt <- function(k12 = 0.2, k21 = 0.2, k10 = 0.1) { curve(Ctissue(x, k12, k21, k10), from=0, to=60, xlab="Time (minutes)", ylab="Drug Concentration (mg/dl)", main="Bi-exponential / Two-compartment Model") curve(CBlood(x, k12, k21, k10), from=0, to=60, add=TRUE, col="darkred") legend("topright", legend=c("Tissue", "Blood"), col=c("black", "red"), lty=1, bty="n") } # bi-exponential model # t_i = [0, 5, 10, 20] # y_ij = Ctissue(t_ij, k12_i, k21_i, k10_i) + e_ij # k_i ~ N_3(mu_k, sigma_k) # e_ij ~ N(0, sigma) # simulate independent predictors independent <- function(n, times=seq(0,20,2)) { return(data.frame( time=rep(times,n), subject=rep(1:n, rep(length(times), n)), sex=rep(c("Female", "Male"), rep(n*length(times)/2,2)))) } # simulate dependant data from a nonlinear mixed effects model dependent <- function(inde, muF = c(0.4, 0.2, 0.2), muM = c(0.2, 0.2, 0.05), sigma_k = diag(c(0.02, 0.02, 0.02)), sigma= 0.005) { inde[c("k12", "k21", "k10")] <- NA np <- length(unique(inde$time)) iF <- which(inde$sex == "Female") nF <- length(unique(inde[iF,]$subject)) kF <- matrix(rnorm(nF*3), nF, 3) kF <- t(t(kF %*% sigma_k) + muF) kF <- kronecker(kF, rep(1,np)) inde[iF, c("k12", "k21", "k10")] <- data.frame(kF) iM <- which(inde$sex == "Male") nM <- length(unique(inde[iM,]$subject)) kM <- matrix(rnorm(nM*3), nM, 3) kM <- t(t(kM %*% sigma_k) + muM) kM <- kronecker(kM, rep(1, np)) inde[iM, c("k12", "k21", "k10")] <- data.frame(kM) inde$y <- Ctissue(inde$time, inde$k12, inde$k21, inde$k10) + rnorm(nrow(inde), 0, sigma) return(inde) } # fit a nonlinear (bi-exponential) mixed effects model fitBiExp <- function(data) { require("nlme") dat <- groupedData(y ~ time | subject, data=data) fit <- nlme(y ~ Ctissue(time, k12, k21, k10), fixed = k12 + k21 + k10 ~ sex, random = pdDiag(k12 + k21 + k10 ~ 1), start = c(0.6, -0.4, 0.3, -0.1, 0.4, -0.35), data=dat) } # plot the data and model fit (including random effects) plotFitBiExp <- function(fit, data) { colFdk <- "#F75E9B" colFlt <- "#FF8ECB" colMdk <- "#577BE6" colMlt <- "#87ABFF" plot(data$time, data$y, col=ifelse(data$sex=="Female", colFlt, colMlt), xlab="Time (minutes)", ylab="Drug Concentration (mg/dl)") pn <- 200 pt <- seq(0, 20, length.out=pn) sub <- unique(data$subject) sex <- data[match(sub, data$subject),]$sex fxf <- fixef(fit) rnf <- ranef(fit) rnf <- rnf[match(as.character(1:length(sub)), rownames(rnf)),] for(i in 1:length(sub)) { sx <- sex[i] k12 <- fxf[1] + rnf[i,1] + if(sx=="Female") 0 else fxf[2] k21 <- fxf[3] + rnf[i,2] + if(sx=="Female") 0 else fxf[4] k10 <- fxf[5] + rnf[i,3] + if(sx=="Female") 0 else fxf[6] curve(Ctissue(x, k12, k21, k10), from=0, to=20, col=if(sx == "Female") colFlt else colMlt, lty=2, add=TRUE) } pM <- predict(fit, data.frame(time=pt, sex="Male", subject=1)) pF <- predict(fit, data.frame(time=pt, sex="Female", subject=1)) lines(pt, pM, col=colMdk, lwd=2) lines(pt, pF, col=colFdk, lwd=2) } @ \normalsize \noindent The figure below illustrates the raw data (points), the fixed effect curves for males (solid blue) and females (solid pink), and the individual random-effect curves for each patient (dashed). \begin{center} <>= data <- dependent(independent(20)) fit <- fitBiExp(data) plotFitBiExp(fit, data) @ \end{center} \end{document}