require("DoseFinding")
## effect curve estimate for linear type models!!
########################################################################
#### Testing function to generate doses and sample size allocs.
genDFdats <- function(model, argsMod, doses, n, sigma, mu = NULL){
nD <- length(doses)
dose <- sort(doses)
if (length(n) == 1) n <- rep(n, nD)
dose <- rep(dose, n)
args <- c(list(dose), argsMod)
mu <- do.call(model, args)
data.frame(dose = dose, resp = mu + rnorm(sum(n), sd = sigma))
}
getDosSampSiz <- function(){
# generate dose levels
mD <- runif(1, 0, 1500)
nD <- max(rpois(1, 5), 4)
p <- rgamma(nD, 3)
p <- cumsum(p/sum(p))
doses <- signif(c(0, mD*p), 3)
# sample size allocations
totSS <- rpois(1, rexp(1, 1/250))
totSS <- max(totSS, 50)
p <- rgamma(nD+1, 3);p <- p/sum(p)
n <- round(p*totSS)
n[n==0] <- rpois(sum(n==0), 1)+1
list(doses=doses, n=n)
}
getDFdataSet <- function(doses, n){
if(missing(doses) & missing(n)){
ll <- getDosSampSiz()
} else {
ll <- list(doses = doses, n=n)
}
e0 <- rnorm(1, 0, 10)
eMax <- rgamma(1, abs(e0)*0.5, 0.5)
sig <- eMax/runif(1, 0.5, 5)
if(runif(1)<0.3){
aa <- genDFdats("betaMod", c(e0 = e0, eMax = eMax, delta1=runif(1, 0.5, 4),
delta2=runif(1, 0.5, 4), scal=1.2*max(ll$doses)),
ll$doses, ll$n, sig)
} else {
aa <- genDFdats("sigEmax", c(e0 = e0, eMax = eMax,
ed50=runif(1, 0.05*max(ll$doses), 1.5*max(ll$doses)),
h=runif(1, 0.5, 4)), ll$doses, ll$n, sig)
}
N <- sum(ll$n)
center <- c("blue", "green", "red", "yellow", "silver")
aa <- data.frame(x= aa$dose, y=aa$resp, center=as.factor(sample(center, N, replace = T)),
age=runif(N, 1, 100))
aa[sample(1:nrow(aa)),]
}
########################################################################
########################################################################
#### Generate data sets and compare results of fitDRModel
#### to the result of nls and lm for AIC function (if these are consistent
#### parameter estimates, residual sum of square and degrees of freedom are
#### consistent) and the vcov function (if these are consistent parameter
#### estimates, RSS, df and gradient are consistent)
########################################################################
########################################################################
#### beta Model
set.seed(2000)
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
bnds <- matrix(c(0.05, 0.05, 6, 6), nrow=2)
fit0 <- fitMod(x, y, datset, model = "betaMod", addCovars = ~1,
addArgs=list(scal=1.2*max(datset$x)), bnds=bnds, start=c(0.6, 0.6))
fitnls <- nls(y~betaMod(x, e0, emax, delta1, delta2, 1.2*max(datset$x)),
start=c(e0=15, emax=14, delta1=0.8, delta2=0.5), data=datset)
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0)
vcov(fitnls)
predict(fit0, predType="effect-curve", se.fit=TRUE)
predict(fit0, predType="full-model", se.fit=TRUE)
TD(fit0, Delta = 1)
# with covariates
fit0 <- fitMod(x, y, datset, model="betaMod", addCovars = ~age+center,
addArgs=list(scal=1.2*max(datset$x)), bnds=bnds)
XX <- model.matrix(~center+age, data=datset)
scl <- 1.2*max(datset$x)
fitnls <- nls(y~cbind(XX, betaMod(x, 0, 1, delta1, delta2, scl)),
data=datset, start=c(delta1=1, delta2=0.2),
algorithm = "plinear")
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", doseSeq = c(0, 100), se.fit=T)
predict(fit0, predType="full-model", se.fit=T,
newdata = data.frame(x = c(0,100), center = as.factor("yellow"), age = 50))
TD(fit0, Delta = 1)
########################################################################
#### emax Model
set.seed(15)
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
bnds <- c(1e-5, max(datset$x))
fit0 <- fitMod(x,y, datset, model="emax", addCovars = ~1, bnds=bnds)
fitnls <- nls(y~emax(x, e0, emax, ed50),
start=c(e0=-1, emax=1.3, ed50=0.1), data=datset)
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", se.fit=T)
predict(fit0, predType="full-model", se.fit=T)
TD(fit0, Delta = 0.005)
# with covariates
fit0 <- fitMod(x,y, datset, model="emax", addCovars = ~age+center, bnds=bnds)
XX <- model.matrix(~center+age, data=datset)
fitnls <- nls(y~cbind(XX, emax(x, 0, 1, ed50)),
data=datset, start=list(ed50=1), algorithm = "plinear")
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", doseSeq = c(0, 100), se.fit=T)
predict(fit0, predType="full-model", se.fit=T,
newdata = data.frame(x = c(0,100), center = as.factor("silver"), age = 50))
TD(fit0, Delta = 0.005)
########################################################################
#### sigEmax Model
## set.seed(25) # example where nls and bndnls find different optimum
set.seed(13)
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
bnds <- matrix(c(1e-5, 1e-5, max(datset$x), 30), nrow=2)
fit0 <- fitMod(x,y, datset, model = "sigEmax", addCovars = ~1, bnds=bnds)
fitnls <- nls(y~sigEmax(x, e0, emax, ed50, h),
start=c(e0=6, emax=17, ed50=240, h=2), data=datset)
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", se.fit=T)
predict(fit0, predType="full-model", se.fit=T)
TD(fit0, Delta = 1)
# with covariates
fit0 <- fitMod(x,y, datset, model="sigEmax", addCovars = ~age+center, bnds=bnds)
XX <- model.matrix(~center+age, data=datset)
fitnls <- nls(y~cbind(XX, sigEmax(x, 0, 1, ed50, h)),
data=datset, start=list(ed50=368, h=2), algorithm = "plinear")
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", doseSeq = c(0, 100), se.fit=T)
predict(fit0, predType="full-model", se.fit=T,
newdata = data.frame(x = c(0,100), center = as.factor("silver"), age = 50))
TD(fit0, Delta = 1)
########################################################################
#### logistic Model
set.seed(200)
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
bnds <- matrix(c(1e-5, 1e-5, max(datset$x), max(datset$x)/2), nrow=2)
fit0 <- fitMod(x,y, datset, model="logistic", addCovars = ~1, bnds=bnds)
fitnls <- nls(y~logistic(x, e0, emax, ed50, delta),
start=c(e0=0, emax=16, ed50=250, delta=90), data=datset)
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", se.fit=T)
predict(fit0, predType="full-model", se.fit=T)
TD(fit0, Delta = 0.5)
# with covariates (example where nls and bndnls find different optima)
fit0 <- fitMod(x,y, datset, model="logistic", addCovars = ~age+center, bnds=bnds)
XX <- model.matrix(~center+age, data=datset)
fitnls <- nls(y~cbind(XX, logistic(x, 0, 1, ed50, delta)),
data=datset, start=list(ed50=220, delta=48), algorithm = "plinear")
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", doseSeq = c(0, 100), se.fit=T)
predict(fit0, predType="full-model", se.fit=T,
newdata = data.frame(x = c(0,100), center = as.factor("silver"), age = 5))
TD(fit0, Delta = 0.02)
########################################################################
#### exponential Model
set.seed(4)
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
bnds <- c(0.1, 2)*max(datset$x)
fit0 <- fitMod(x,y, datset, model = "exponential", addCovars = ~1, bnds=bnds)
fitnls <- nls(y~exponential(x, e0, e1, delta),
start=coef(fit0), data=datset)
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", se.fit=T)
predict(fit0, predType="full-model", se.fit=T)
TD(fit0, Delta = 0.1)
# with covariates
bnds <- c(0.1, 2)*max(datset$x)
fit0 <- fitMod(x,y, datset, model = "exponential", addCovars = ~age+center,
bnds=bnds)
XX <- model.matrix(~center+age, data=datset)
fitnls <- nls(y~cbind(XX, exponential(x, 0, 1, delta)),
data=datset, start=c(delta=450), algorithm = "plinear")
AIC(fit0)
AIC(fitnls)
summary(fit0)
summary(fitnls)
vcov(fit0 )
vcov(fitnls)
predict(fit0, predType="effect-curve", doseSeq = c(0, 100), se.fit=T)
predict(fit0, predType="full-model", se.fit=T,
newdata = data.frame(x = c(0,100), center = as.factor("blue"), age = 50))
TD(fit0, Delta = 0.1)
########################################################################
#### linear model
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
fit0 <- fitMod(x,y, datset, model = "linear", addCovars = ~1)
fitlm <- lm(y~x, data=datset)
AIC(fit0)
AIC(fitlm)
summary(fit0)
summary(fitlm)
vcov(fit0 )
vcov(fitlm)
predict(fit0, predType="effect-curve", se.fit=T)
TD(fit0, Delta = 1)
# with covariates
fit0 <- fitMod(x,y, datset, model = "linear", addCovars = ~age+center)
fitlm <- lm(y~x+age+center, data=datset)
AIC(fit0)
AIC(fitlm)
summary(fit0)
summary(fitlm)
vcov(fit0 )
vcov(fitlm)
predict(fit0, predType="effect-curve", se.fit=T)
predict(fit0, predType = "f", se.fit = T,
newdata = data.frame(x=c(0,1,2,100), age = 30, center = as.factor("blue")))
predict(fitlm, se.fit = T,
newdata = data.frame(x=c(0,1,2,100), age = 30, center = as.factor("blue")))
TD(fit0, Delta = 1)
########################################################################
#### linlog model
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
off <- 0.05*max(datset$x)
# without covariates
fit0 <- fitMod(x,y, datset, model = "linlog", addCovars = ~1,addArgs=list(off=off))
fitlm <- lm(y~log(x+off), data=datset)
AIC(fit0)
AIC(fitlm)
summary(fit0)
summary(fitlm)
vcov(fit0 )
vcov(fitlm)
predict(fit0, predType="effect-curve", se.fit=T) ## bug ##
TD(fit0, Delta = 1)
# with covariates
fit0 <- fitMod(x,y, datset, model = "linlog", addCovars = ~age+center,
addArgs=list(off=off))
fitlm <- lm(y~log(x+off)+age+center, data=datset)
AIC(fit0)
AIC(fitlm)
summary(fit0)
summary(fitlm)
vcov(fit0 )
vcov(fitlm)
predict(fit0, predType = "f", se.fit = T, ## degrees of freedom wrong ##
newdata = data.frame(x=c(0,1,2,100), age = 35, center = as.factor("blue")))
predict(fitlm, se.fit = T,
newdata = data.frame(x=c(0,1,2,100), age = 35, center = as.factor("blue")))
TD(fit0, Delta = 1)
########################################################################
#### quadratic model
ll <- getDosSampSiz()
datset <- getDFdataSet(ll$doses, ll$n)
# without covariates
fit0 <- fitMod(x,y, datset, model = "quadratic", addCovars = ~1)
fitlm <- lm(y~x+I(x^2), data=datset)
AIC(fit0)
AIC(fitlm)
summary(fit0)
summary(fitlm)
vcov(fit0 )
vcov(fitlm)
predict(fit0, predType="effect-curve", se.fit=T)
predict(fit0, predType="full-model", se.fit=T,
newdata=data.frame(x=c(0, 10, 100)))
predict(fitlm, se.fit=T,
newdata=data.frame(x=c(0, 10, 100)))
TD(fit0, Delta = 1)
# with covariates
fit0 <- fitMod(x,y, datset, model = "quadratic", addCovars = ~age+center)
fitlm <- lm(y~x+I(x^2)+age+center, data=datset)
AIC(fit0)
AIC(fitlm)
summary(fit0)
summary(fitlm)
vcov(fit0 )
vcov(fitlm)
predict(fit0, predType = "f", se.fit = T,
newdata=data.frame(x=c(0, 10, 100), age = 30, center = as.factor("blue")))
predict(fitlm, se.fit = T,
newdata=data.frame(x=c(0, 10, 100), age = 30, center = as.factor("blue")))
TD(fit0, Delta = 0.1)
########################################################################
## ensure that predict with no argument uses the original data not the
## sorted data that were used for fitting
data(IBScovars)
ff <- fitMod(dose, resp, data=IBScovars, model="quadratic",
addCovars = ~gender)
## should be all zero
predict(ff, predType = "ls-means")-
predict(ff, predType = "ls-means", doseSeq = IBScovars[,3])
predict(ff, predType = "full-model")-
predict(ff, predType = "full-model", newdata = IBScovars[,-2])
predict(ff, predType = "effect-curve")-
predict(ff, predType = "effect-curve", doseSeq = IBScovars[,3])
ff2 <- fitMod(dose, resp, data=IBScovars, model="quadratic")
## should be all zero
predict(ff2, predType = "ls-means")-
predict(ff2, predType = "ls-means", doseSeq = IBScovars[,3])
predict(ff2, predType = "full-model")-
predict(ff2, predType = "full-model", newdata = IBScovars[,-2])
predict(ff2, predType = "effect-curve")-
predict(ff2, predType = "effect-curve", doseSeq = IBScovars[,3])
dose <- unique(IBScovars$dose)
ord <- c(2,4,1,3,5)
mns <- tapply(IBScovars$resp, IBScovars$dose, mean)[ord]
ff3 <- fitMod(dose, mns, S=diag(5), model="quadratic", type = "general")
predict(ff3, predType = "ls-means")-
predict(ff3, predType = "ls-means", doseSeq = dose)
predict(ff3, predType = "effect-curve")-
predict(ff3, predType = "effect-curve", doseSeq = dose)
########################################################################
## ensure that S is also sorted when the dose is not entered sorted
dose <- sort(unique(IBScovars$dose))
mns <- tapply(IBScovars$resp, IBScovars$dose, mean)
S <- c(1000,1,1,1,1)*diag(5)
ff1 <- fitMod(dose, mns, S = S, model="linear", type="general")
## fit unsorted
dose <- unique(IBScovars$dose)
ord <- c(2,4,1,3,5)
mns <- tapply(IBScovars$resp, IBScovars$dose, mean)[ord]
ff2 <- fitMod(dose, mns, S = S, model="linear", type="general")
ff3 <- fitMod(dose, mns, S = S[ord,ord], model="linear", type="general")
## coef(ff1) & coef(ff3) should be equal
coef(ff1)
coef(ff3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.