Dear HS !
thank you for your input and the references!
below some code to play around; it is an optimal contrast problem – the power to detect a deviation from dose-proportionality depend on the unknown monotone dose-response shape. Theoretically, the hypothesis regarding dose-proportionality should be specified in detail to “maximise power”. For example, the hypothesis
H0: all 4 doses are within the proportional range
H1: 3 highest doses are outside the proportional range
would suggest application of the reverse helmert contrast. Additionally, one-sided step-down testing allows a formal estimation of the minimum dose where we can not reject the null hypothesis of dose-proportionality – similar as estimation of the minimum effective dose.
best regards
martin
rm(list=ls(all=TRUE))
# function to simulate power for an overall assessment of dose-proportionality based on the power-law model and 4 dose levels
sim.dose4prop <- function(nsim, cv, n, a, b, mindose, maxdose, range, alpha=0.05){
# input parameters:
# nsim ... number of simulation runs
# cv ... coefficient of variation of random variables (e.g. AUCs)
# n ... n per dose level
# a, b ... parameters for power law model: AUC = a*dose^b
# mindose ... lowest dose
# maxdose ... highest dose
# range ... doses (index) for which parameter b of power law model is applied (e.g. b=1.1 and range=c(2:4) specifies that doses 2 to 4 are outside proportional range)
# alpha ... type I error, two-sided (default = 0.05)
# function to define AUC using power law model
powermod <- function(dose, a, b, range){return(c(a*dose[-range]^1, a*dose[range]^b))}
# function to generate log-normal distributed RVs (e.g. AUCs) with identical CVs per dose. See e.g. Jaki et al. (2009) for more detail.
# Jaki T, Wolfsegger MJ, Ploner M (2009). Confidence intervals for ratios of AUCs in the case of serial sampling: A comparison of seven methods.
# Pharmaceutical Statistics, 8(1):12-24. DOI:10.1002/pst.321.
data.gen <- function(dose, mue, n, cv){
k <- length(dose)
sd.log <- sqrt(log(1+cv^2))
res <- as.data.frame(matrix(nrow=k*n, ncol=6))
res[,1] <- sort(rep(dose,n))
res[,2] <- sort(rep(mue,n))
res[,3] <- log(res[,2]) - 1/2 * log(1+cv^2)
res[,4] <- sd.log
res[,5] <- rlnorm(n=k*n, meanlog=res[,3], sdlog=res[,4])
res[,6] <- res[,5]/res[,1]
colnames(res) <- c('dose', 'mue', 'mue.log', 'sd.log', 'auc', 'auc.adj')
return(res)
}
# definition of contrasts see e.g. Bretz and Hothorn (2003) for more detail
# Bretz F and Hothorn LA (2003). Statistical Analysis of Monotone or Non-monotone Dose–Response Data from In Vitro Toxicological Assays.
# ATLA-ALTERN LAB ANIM 31:81-96 Suppl. 1 JUN 2003. URL: http://ecvam.jrc.it/publication/Bretz-Hothorn.pdf last visited at 2009-12-27
numdose <- 4
contr.lowhigh <- c(-1,0,0,1) # t-test of highest versus lowest dose (pairwise contrast)
contr.helmert <- c(-1,-1,-1,3) # helmert contrast
contr.revhelm <- c(-3,1,1,1) # reverse helmert contrast
contr.linearc <- c(-3,-1,1,3) # linear contrast
contr.jonckht <- c(-1, -0.67, -0.17, 1.83) # Jonkheere-Terpstra analogon contrast
contr.maximum <- c(-0.87, -0.13, 0.13, 0.87) # Maximum contrast
contr.linr24t <- c(-12, -2, 2, 12) # Linear 2–4 contrast
# define doses: to be equidistant on log-scale
step <- (log(maxdose)-log(mindose))/(numdose-1)
log.dose <- c(log(mindose), rep(NA, numdose-1))
for(i in 2:numdose){log.dose[i] <- log.dose[i-1]+step}
dose <- exp(log.dose)
mue <- powermod(dose=dose, a=a, b=b, range=range)
k <- length(dose)
power <- rep(0,9)
for(i in 1:nsim){
# generate log-normal distributed random variables
data <- data.gen(dose=dose, mue=mue, n=n, cv=cv)
# log-log regression approach: two-sided CI
res1 <- confint(lm(log(data$auc) ~ log(data$dose)), level=1-alpha)[2,]
# ANOVA approach
temp <- summary(lm(log(data$auc.adj) ~ as.factor(data$dose)))
res2 <- 1-pf(q=temp$fstatistic[1], df1=temp$fstatistic[2], df2=temp$fstatistic[3])
# contrast approaches: two-sided hypothesis
xq <- tapply(log(data$auc.adj), data$dose, mean)
sd.pooled <- sqrt(sum((n-1)*tapply(log(data$auc.adj), data$dose, var))/(n*k-k))
df <- n*k - k
nom3 <- sum(contr.lowhigh*xq)
den3 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.lowhigh)^2))
nom4 <- sum(contr.helmert*xq)
den4 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.helmert)^2))
nom5 <- sum(contr.revhelm*xq)
den5 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.revhelm)^2))
nom6 <- sum(contr.linearc*xq)
den6 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.linearc)^2))
nom7 <- sum(contr.jonckht*xq)
den7 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.jonckht)^2))
nom8 <- sum(contr.maximum*xq)
den8 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.maximum)^2))
nom9 <- sum(contr.linr24t*xq)
den9 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.linr24t)^2))
res3 <- 2*(1-pt(q=abs(nom3/den3), df=df))
res4 <- 2*(1-pt(q=abs(nom4/den4), df=df))
res5 <- 2*(1-pt(q=abs(nom5/den5), df=df))
res6 <- 2*(1-pt(q=abs(nom6/den6), df=df))
res7 <- 2*(1-pt(q=abs(nom7/den7), df=df))
res8 <- 2*(1-pt(q=abs(nom8/den8), df=df))
res9 <- 2*(1-pt(q=abs(nom9/den9), df=df))
if(res1[1] > 1 | res1[2] < 1){power[1] <- power[1] + 1}
if(res2 < alpha){power[2] <- power[2] + 1}
if(res3 < alpha){power[3] <- power[3] + 1}
if(res4 < alpha){power[4] <- power[4] + 1}
if(res5 < alpha){power[5] <- power[5] + 1}
if(res6 < alpha){power[6] <- power[6] + 1}
if(res7 < alpha){power[7] <- power[7] + 1}
if(res8 < alpha){power[8] <- power[8] + 1}
if(res9 < alpha){power[9] <- power[9] + 1}
}
res <- rep(NA, 11)
res[1] <- mindose
res[2] <- maxdose
res[c(3:11)] <- power/nsim
names(res) <- c('mindose', 'maxdose', 'log-log reg', 'anova', 't-test', 'helmert', 'revhelm', 'linearc', 'jonckht', 'maximum', 'linr24t')
return(res)
}
# set parameters for simulations
nsim <- 1E4
cv <- 0.5
n <- 10
a <- 1
# dose-proportionality for complete dose range (simulate type I error)
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1, mindose=10, maxdose=20, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1, mindose=10, maxdose=100, range=c(1:4))
# all 4 doses outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(1:4))
# last 3 doses outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(2:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(2:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(2:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(2:4))
# last 2 doses outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(3:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(3:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(3:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(3:4))
# only highest dose outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(4))
|