## (i)
## example from Biometrics paper p. 743
data(biom)
models <- list(linear = NULL, linlog = NULL, emax = 0.2,
exponential = c(0.279,0.15), quadratic = c(-0.854,-1))
dfe <- MCPMod(resp ~ dose, biom, models, alpha = 0.05, doseEstPar = 0.05,
pVal = TRUE, selModel = "maxT", doseEst = "MED2old",
clinRel = 0.4, off = 1)
## detailed information is available via summary
summary(dfe)
## plots data with selected model function
plot(dfe)
## example with IBS data
## now explicitly calculate critical value and perform
## model averaging based on the AIC
data(IBScovars)
models <- list(emax = 0.2, quadratic = -0.2, linlog = NULL)
dfe2 <- MCPMod(resp ~ dose, IBScovars, models, alpha = 0.05, critV = TRUE,
pVal = TRUE, selModel = "aveAIC",
clinRel = 0.25, off = 1)
dfe2
## illustrate some methods for MCPMod objects
summary(dfe2)
plot(dfe2, complData = TRUE)
plot(dfe2, CI = TRUE, clinRel = TRUE)
## use different optimizer (and include some models
## not converging using nls)
models <- list(quadratic = -0.2, linlog = NULL, betaMod = c(1,1),
sigEmax = rbind(c(0.2,1), c(1, 3)))
dfe3 <- MCPMod(resp ~ dose, IBScovars, models, alpha = 0.05, pVal = TRUE,
selModel = "aveAIC", clinRel = 0.25, off = 1,
scal = 6, optimizer = "bndnls")
plot(dfe3)
## use additional linear covariates
data(IBScovars)
models <- list(emax = 0.2, quadratic = -0.2, linlog = NULL)
dfe4 <- MCPMod(resp ~ dose, IBScovars, models, addCovars = ~gender,
alpha = 0.05, pVal = TRUE,
selModel = "aveAIC", clinRel = 0.25, off = 1)
plot(dfe4, CI = TRUE) # plot method now only plots the effect curves
## simulate dose-response data
set.seed(1)
dfData <- genDFdata(model = "emax", argsMod = c(e0 = 0.2, eMax = 1,
ed50 = 0.05), doses = c(0,0.05,0.2,0.6,1), n=20, sigma=0.5)
models <- list(emax = 0.1, logistic = c(0.2, 0.08),
betaMod = c(1, 1))
## hand over starting values manually
start <- list(emax = c(ed50 = 0.1), logistic = c(ed50=0.05, delta =
0.1), betaMod = c(delta1 = 0.5, delta2 = 0.5))
dfe5 <- MCPMod(resp ~ dose, dfData, models, clinRel = 0.4, critV = 1.891,
scal = 1.5, start = start)
## (ii) Example for constructing a model list
## Contrasts to be included:
## Model guesstimate(s) for stand. model parameter(s) (name)
## linear -
## linear in log -
## Emax 0.2 (ED50)
## Emax 0.3 (ED50)
## exponential 0.7 (delta)
## quadratic -0.85 (delta)
## logistic 0.4 0.09 (ED50, delta)
## logistic 0.3 0.1 (ED50, delta)
## betaMod 0.3 1.3 (delta1, delta2)
## sigmoid Emax 0.5 2 (ED50, h)
##
## For each model class exactly one list entry needs to be used.
## The names for the list entries need to be written exactly
## as the model functions ("linear", "linlog", "quadratic", "emax",
## "exponential", "logistic", "betaMod", "sigEmax").
## For models with no parameter in the standardized model just NULL is
## specified as list entry.
## For models with one parameter a vector needs to be used with length
## equal to the number of contrasts to be used for this model class.
## For the models with two parameters in the standardized model a vector
## is used to hand over the contrast, if it is desired to use just one
## contrast. Otherwise a matrix with the guesstimates in the rows needs to
## be used. For the above example the models list has to look like this
list(linear = NULL, linlog = NULL, emax = c(0.2, 0.3),
exponential = 0.7, quadratic = -0.85, logistic =
rbind(c(0.4, 0.09), c(0.3, 0.1)),
betaMod = c(0.3, 1.3), sigEmax = c(0.5, 2))
## Additional parameters (which will not be estimated) are the "off"
## parameter for the linlog model and the "scal" parameter for the
## beta model, which are not handed over in the model list.
## (iii) example for incorporation of a usermodel
## simulate dose response data
dats <- genDFdata("sigEmax", c(e0 = 0, eMax = 1, ed50 = 2, h = 2),
n = 50, sigma = 1, doses = 0:10)
## define usermodel
userMod <- function(dose, a, b, d){
a + b*dose/(dose + d)
}
## define gradient
userModGrad <-
function(dose, a, b, d) cbind(1, dose/(dose+d), -b*dose/(dose+d)^2)
## name starting values for nls
start <- list(userMod=c(a=1, b=2, d=2))
models <- list(userMod=c(0, 1, 1), linear = NULL)
dfe6 <- MCPMod(resp ~ dose, dats, models, clinRel = 0.4, selModel="AIC",
start = start, uGrad = userModGrad)
## (iv) Contrast matrix and critical value handed over and not calculated
## simulate dose response data
dat <- genDFdata(mu = (0:4)/4, n = 20,
sigma = 1, doses = (0:4)/4)
## construct optimal contrasts and critical value with planMM
doses <- (0:4)/4
mods <- list(linear = NULL, quadratic = -0.7)
pM <- planMM(mods, doses, 20)
fit <- MCPMod(resp ~ dose, dat, models = NULL, clinRel = 0.3,
contMat = pM$contMat, critV = pM$critVal)
## (v) Use self constructed contrasts (and hand them over)
mu1 <- c(1, 2, 2, 2, 2)
mu2 <- c(1, 1, 2, 2, 2)
mu3 <- c(1, 1, 1, 2, 2)
mMat <- cbind(mu1, mu2, mu3)
dimnames(mMat)[[1]] <- doses
pM <- planMM(muMat = mMat, doses = doses, n = 20, cV = FALSE)
## use MCPMod only for testing part (see also MCPtest function)
fit <- MCPMod(resp ~ dose, dat, models = NULL, contMat = pM$contMat,
pVal = TRUE, testOnly = TRUE, critV=TRUE)
summary(fit)Run the code above in your browser using DataLab