Nothing
#' Prepare General Data for the MCPMod Function
#'
#' This function serves as an alternative for using the MCPModGen function
#' directly for general data. The function returns the estimates for \eqn{\mu} and
#' \eqn{S}, which are needed for MCPMod.
#'
#' @param family A character string containing the error distribution to be used
#' in the model.
#' @param link A character string for the model link function.
#' @param w Either a numeric vector of the same length as dose and resp, or a
#' character vector denoting the column name in the data.
#' @param offset Either a numeric vector of the same length as dose and resp, or a
#' character vector denoting the column name in the data.
#' @param ... Additional arguments to be passed to \code{glm} or \code{glm.nb}.
#' This is especially useful when a fitting error is returned. In these cases,
#' it may be useful to supply a \code{start} vector for the parameters.
#' @inheritParams MCPModGen
#' @return A list containing the \eqn{\mu} vector and \eqn{S} matrix.
#' @examples
#' # Analyze the binary migraine data from the DoseFinding package.
#' data(migraine)
#' models = Mods(linear = NULL, emax = 1, quadratic = c(-0.004), doses = migraine$dose)
#'
#' # Now analyze using binomial weights
#' PFrate <- migraine$painfree/migraine$ntrt
#' migraine$pfrat = migraine$painfree / migraine$ntrt
#' muS = prepareGen("binomial", "logit", w = "ntrt", dose = "dose",
#' resp = "pfrat", data = migraine)
#' ## Look at the elements of muS
#' muS
#' MCPMod(muS$data$dose, muS$data$resp, models = models, S = muS$S,
#' type = "general", selModel = "aveAIC",Delta = 0.2)
#' @export
prepareGen <- function(family = c("negative binomial", "binomial", "poisson"),
link = c("log", "logit", "probit", "cauchit", "cloglog",
"identity", "log risk ratio", "risk ratio",
"sqrt"), w = NULL, dose, resp, data = NULL,
addCovars = ~1, placAdj = FALSE,
offset = NULL, ...) {
## Format user data
if (!is.null(data)) {
if (!inherits(data, "data.frame")) {
stop("data must be of class \"data.frame\"")
}
if (!inherits(dose, "character") | !inherits(resp, "character")) {
stop("dose and resp must be character vectors")
}
dose.vec <- data[[dose]]
resp.vec <- data[[resp]]
if (!is.null(w)) {
if (!inherits(w, "character")) {
stop("w must be a character vector")
}
w.vec <- data[[w]]
}
} else {
if (length(dose) != length(resp)) {
stop("dose and resp must be of equal length")
}
dose.vec <- dose
resp.vec <- resp
if (!is.null(w)) {
if (length(w) != length(resp)) {
stop("w must be of equal length as dose and resp")
}
w.vec <- w
}
if (addCovars != ~1) {
stop("must supply data when addCovars is given")
}
}
if (addCovars == ~1) {
if (!is.null(w)) {
dat.glm <- data.frame(dose = dose.vec, resp = resp.vec, w = w.vec)
} else {
dat.glm <- data.frame(dose = dose.vec, resp = resp.vec)
}
} else {
dat.glm <- data
## Need to rename dose and resp columns
which.dose <- which(names(dat.glm) == dose)
which.resp <- which(names(dat.glm) == resp)
names(dat.glm)[which.dose] = "dose"
names(dat.glm)[which.resp] = "resp"
}
if (is.null(offset)) {
## Do nothing
} else if (!is.null(data) & inherits(offset, "character")) {
offset <- data[[offset]]
} else if (length(offset) != nrow(dat.glm)) {
stop("offset must be the same length as data")
}
doses <- sort(unique(dat.glm$dose))
n.doses <- length(doses)
family <- match.arg(family)
link.in <- match.arg(link)
## Negative binomial, binomial, and Poisson are similar, so comments are only
## included for the negative binomial chunk
if (family == "negative binomial") {
if (!is.null(w)) {
message("ignoring w argument for negative binomial")
}
## For log risk ratio and risk ratio, we use a log link We create the formula then
## call MASS::glm.nb. Log risk ratio is identical to placebo-adjusted
## on log-scale, so fit that addCovars[2] gets the arguments part of the addCovars
if (link.in == "log risk ratio" | link.in == "risk ratio") {
if (is.null(offset)) {
## If there is no offset
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
} else {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2],
"+ offset(log(offset))", sep = ""))
}
# call.glm <- list(form, dat.glm, "log")
# names(call.glm) <- c("formula", "data", "link")
# glm.fit <- do.call(MASS::glm.nb, call.glm)
glm.fit = MASS::glm.nb(form, data = dat.glm, link = "log", ...)
mu.hat <- coef(glm.fit)
S.hat <- vcov(glm.fit)
if (link == "risk ratio") {
## Convert to risk-ratio using Delta method
mu.hat <- exp(mu.hat)
S.hat <- diag(mu.hat) %*% S.hat %*% diag(mu.hat)
}
mu.hat <- mu.hat[1:n.doses]
S.hat <- S.hat[1:n.doses, 1:n.doses]
} else {
if (placAdj == FALSE) {
if (is.null(offset)) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], "-1",
sep = ""))
} else {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2],
"-1 + offset(log(offset))", sep = ""))
}
call.glm <- list(form, dat.glm, link.in)
names(call.glm) <- c("formula", "data", "link")
glm.fit <- do.call(MASS::glm.nb, call.glm)
## I would like to use the MASS::glm.nb version so I can pass the ...
## argument, or find a way to pass ... through do.call. However,
## MASS::glm.nb fails to pass link.in inside the Poisson call
# glm.fit <- MASS::glm.nb(form, data = dat.glm, link = link.in, ...)
mu.hat <- coef(glm.fit)[1:n.doses]
S.hat <- vcov(glm.fit)[1:n.doses, 1:n.doses]
} else if (placAdj) {
if (is.null(offset)) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
} else {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2],
"+ offset(log(offset))", sep = ""))
}
call.glm <- list(form, dat.glm, link.in)
names(call.glm) <- c("formula", "data", "link")
glm.fit <- do.call(MASS::glm.nb, call.glm)
# glm.fit <- MASS::glm.nb(form, data = dat.glm, link = link.in, ...)
mu.hat <- coef(glm.fit)[2:n.doses]
S.hat <- vcov(glm.fit)[2:n.doses, 2:n.doses]
doses <- doses[2:n.doses]
} else {
stop("placAdj must be logical")
}
}
} else if (family == "binomial") {
if (link.in == "log risk ratio" | link.in == "risk ratio") {
if (is.null(w)) {
risks.vec = rep(NA, n.doses)
k <- 1
for (i in doses) {
risks.vec[k] <- mean(dat.glm$resp[dat.glm$dose == i])
k <- k + 1
}
emp.rr <- risks.vec/risks.vec[1]
if (length(labels(terms(addCovars))) > 0) {
emp.rr <- c(emp.rr, rep(1, length(labels(terms(addCovars)))))
}
## Calculate log-risk ratio glm.fit <- glm(resp ~ factor(dose), data = dat.glm,
## family = binomial(link = 'log'), start = log(c(mean(dat.glm$resp[dat.glm$dose
## == 0]), emp.rr[-1])))
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
# call.glm <- list(form, dat.glm, binomial(link = "log"),
# log(c(mean(dat.glm$resp[dat.glm$dose == 0]), emp.rr[-1])))
# names(call.glm) <- c("formula", "data", "family", "start")
# glm.fit <- do.call(glm, call.glm)
glm.fit <- glm(form, data = dat.glm, family = binomial(link = "log"),
start = log(c(mean(dat.glm$resp[dat.glm$dose == 0]),
emp.rr[-1])))
mu.hat <- coef(glm.fit)
S.hat <- vcov(glm.fit)
if (link.in == "risk ratio") {
## Convert to risk-ratio using Delta method
mu.hat <- exp(mu.hat)
S.hat <- diag(mu.hat) %*% S.hat %*% t(diag(mu.hat))
}
mu.hat <- mu.hat[1:n.doses]
S.hat <- S.hat[1:n.doses, 1:n.doses]
} else {
risks.vec <- rep(NA, n.doses)
k <- 1
for (i in doses) {
risks.vec[k] <- dat.glm$resp[dat.glm$dose == i]
k <- k + 1
## there is a faster way to do this, but this ensures the order otherwise use
## emp.rr = dat.glm$resp/dat.glm$resp[dat.glm$dose == i]
}
emp.rr <- risks.vec/risks.vec[1]
if (length(labels(terms(addCovars))) > 0) {
emp.rr <- c(emp.rr, rep(1, length(labels(terms(addCovars)))))
}
# glm.fit <- glm(resp ~ factor(dose), data = dat.glm, family = binomial(link =
# 'log'), weights = w, start = log(c(dat.glm$resp[dat.glm$dose == 0],
# emp.rr[-1])))
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
# call.glm <- list(form, dat.glm, binomial(link = "log"),
# log(c(dat.glm$resp[dat.glm$dose == 0], emp.rr[-1])), dat.glm$w)
# names(call.glm) <- c("formula", "data", "family", "start", "weights")
# glm.fit <- do.call(glm, call.glm)
glm.fit <- glm(form, data = dat.glm, family = binomial(link = "log"),
weights = dat.glm$w,
start = log(c(dat.glm$resp[dat.glm$dose == 0], emp.rr[-1])))
mu.hat <- coef(glm.fit)
S.hat <- vcov(glm.fit)
if (link.in == "risk ratio") {
## Convert to risk-ratio using Delta method
mu.hat <- exp(mu.hat)
S.hat <- diag(mu.hat) %*% S.hat %*% t(diag(mu.hat))
}
mu.hat <- mu.hat[1:n.doses]
S.hat <- S.hat[1:n.doses, 1:n.doses]
}
} else {
if (placAdj == FALSE) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], "-1",
sep = ""))
glm.fit <- glm(form, data = dat.glm, weights = w,
family = binomial(link = link.in), ...)
mu.hat <- coef(glm.fit)[1:n.doses]
S.hat <- vcov(glm.fit)[1:n.doses, 1:n.doses]
} else if (placAdj) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
glm.fit <- glm(form, data = dat.glm, weights = w,
family = binomial(link = link.in), ...)
mu.hat <- coef(glm.fit)[2:n.doses]
S.hat <- vcov(glm.fit)[2:n.doses, 2:n.doses]
doses <- doses[2:n.doses]
} else {
stop("placAdj must be logical")
}
}
} else if (family == "poisson") {
if (!is.null(w)) {
message("ignoring w argument for Poisson")
}
if (link.in == "log risk ratio" | link.in == "risk ratio") {
if (is.null(offset)) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
} else {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2],
" + offset(log(offset))", sep = ""))
}
glm.fit <- glm(form, data = dat.glm, family = poisson(link = "log"),
...)
mu.hat <- coef(glm.fit)
S.hat <- vcov(glm.fit)
if (link.in == "risk ratio") {
## Convert to risk-ratio using Delta method
mu.hat <- exp(mu.hat)
S.hat <- diag(mu.hat) %*% S.hat %*% diag(mu.hat)
}
mu.hat <- mu.hat[1:n.doses]
S.hat <- S.hat[1:n.doses, 1:n.doses]
} else {
if (placAdj == FALSE) {
if (is.null(offset)) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], "-1",
sep = ""))
} else {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2],
"-1 + offset(log(offset))", sep = ""))
}
glm.fit <- glm(form, data = dat.glm, family = poisson(link = link.in),
...)
mu.hat <- coef(glm.fit)[1:n.doses]
S.hat <- vcov(glm.fit)[1:n.doses, 1:n.doses]
} else if (placAdj) {
form <- as.formula(paste("resp~factor(dose)+", addCovars[2], sep = ""))
glm.fit <- glm(form, data = dat.glm, family = poisson(link = link.in),
...)
mu.hat <- coef(glm.fit)[2:n.doses]
S.hat <- vcov(glm.fit)[2:n.doses, 2:n.doses]
doses <- doses[2:n.doses]
} else {
stop("placAdj must be logical")
}
}
} else {
stop("invalid family argument")
}
data.df <- data.frame(dose = doses, resp = mu.hat)
return.list <- list(data = data.df, S = S.hat)
return(return.list)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.