R/aftfit.R
In jaredhuling/aftiv: Instrumental Variable Estimation for Time-to-Event Data
Defines functions
aftfit
repFitAFT
fitAFT
bootstrap.model
Documented in
aftfit
#' Instrumental variable estimation via semiparametric rank-based AFT models
#'
#' @description Fitting function for instrumental variable estimation under the semiparametric AFT model
#' @param formula a formula object, with the response on the left of a tilde operator, and the covariates.
#' The response must be a survival object as returned by the \code{Surv()} function.
#' @param data a \code{data.frame} which contains the variables named in the \code{formula}.
#' @param instrument a vector of length equal to the number of rows in \code{data} containing the instrumental variable
#' @param confounded.x.names name of the exposure variable (endogenous variable) of interest
#' @param method one of:
#' \itemize{
#' \item{\code{"AFT-IPCW"}}{ The proposed IV method that uses inverse probability of censoring weighting }
#' \item{\code{"AFT"}} { Plain rank-based AFT method that does not account for unmeasured confounding}
#' \item{\code{"AFT-IV"}} { biased AFT method that acccounts for unmeasured confounding, but not the non-ignorable censoring
#' (and is hence biased) }
#' \item{\code{"AFT-2SLS"}} { Two-stage AFT method which accounts for unmeasured confounding but requires correct specification
#' of relationship between IV and exposure }
#' }
#' @param smoothed should a smoothed version of the estimating equation be used? Experimental.
#' @param weights observation weights. Not used currently!
#' @param bootstrap logical variable indicating whether or not to run bootstrap
#' @param boot.method bootstrap method to use. \code{"ls"} uses the LS method of Zeng and Lin (2008),
#' \code{"sv"} uses the SV method of Zeng and Lin (2008). \code{"full.bootstrap"} uses a full boostrap where resampled estimating equations
#' are solved for each iteration (this method is very slow!)
#' @param B number of bootstrap iterations
#' @param dependent.censoring for the \code{"AFT-IPCW"} method, should the censoring model allow for dependence on covariates? If \code{TRUE}
#' a Cox model will be fit for the censoring distribution
#' @param na.action how should missing data be treated? See help file for \code{\link[stats]{lm}} for more details
#' @param return.data should the data be returned?
#' @param init vector equal to the length of the number of parameters in the model used to initialize estimation
#' @param tol positive tolerance threshold. Smaller values indicate more precise solutions are required
#' @param maxit maximum number of restarts for \code{BB} algorithm
#' @param verbose should messages be printed? logical variable
#' @param BB.control control list for \code{BBsolve()} function. See \code{\link[BB]{BBsolve}} for more information
#' @param ... not used
#' @import ggplot2
#' @import stats
#' @export
#'
#'
#' @examples
#' library(aftiv)
#'
#' set.seed(1)
#' true.beta <- c(1,-0.5,-0.5,0.25,0,0,-0.75,0.75)
#' dat <- simIVMultivarSurvivalData(500,1,1,-1,1,true.beta,num.confounded = 1,
#' confounding.function = "exp")
#'
#' df <- data.frame(dat$survival[c("delta", "log.t")], dat$X)
#' Z <- dat$Z
#'
#' system.time(aftf <- aftfit(Surv(log.t, delta) ~ ., data = df,
#' instrument = Z,
#' confounded.x.names = "X1",
#' method = c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW"),
#' boot.method = "ls", verbose = 1,
#' B = 200L, smoothed = FALSE,
#' bootstrap = TRUE))
#'
#' summary(aftf)
aftfit <- function(formula,
data,
instrument,
confounded.x.names,
method = c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW"),
smoothed = FALSE,
weights = NULL,
bootstrap = FALSE,
boot.method = c("ls", "sv", "full.bootstrap"),
B = 1000L,
dependent.censoring = FALSE,
na.action,
init = NULL,
return.data = FALSE,
tol = 1e-5,
maxit = 10L,
verbose = 0,
BB.control = NULL,
...) {
method <- match.arg(method, several.ok = TRUE)
boot.method <- match.arg(boot.method)
Call <- match.call()
# create a call to model.frame() that contains the formula (required)
# and any other of the relevant optional arguments
# then evaluate it in the proper frame
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
names(Call), nomatch=0)
if (indx[1] == 0) stop("A formula argument is required")
temp <- Call[c(1,indx)] # only keep the arguments we wanted
temp[[1]] <- as.name('model.frame') # change the function called
temp$formula <- if(missing(data)) terms(formula)
else terms(formula, data=data)
mf <- eval(temp, parent.frame())
if (nrow(mf) == 0) stop("No (non-missing) observations")
Terms <- terms(mf)
Y <- model.extract(mf, "response")
if (!inherits(Y, "Surv")) stop("Response must be a survival object")
type <- attr(Y, "type")
if (type!='right' && type!='counting')
stop(paste("AFT model doesn't support \"", type,
"\" survival data", sep=''))
weights <- model.weights(mf)
nobs <- nrow(Y) #remember this before any time transforms
types <- c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW")
funcs <- c("AFTScorePre", "AFTivScorePre", "AFT2SLSScorePre", "AFTivIPCWScorePre")
funcs.sm <- c("AFTScoreSmoothPre", "AFTivScoreSmoothPre", "AFT2SLSScoreSmoothPre", "AFTivIPCWScoreSmooth")
contrast.arg <- NULL #due to shared code with model.matrix.coxph
X <- model.matrix(Terms, mf, contrasts=contrast.arg)
#print(head(X))
# drop the intercept after the fact, and also drop strata if necessary
Xatt <- attributes(X)
adrop <- 0
xdrop <- Xatt$assign %in% adrop
#xdrop <- Xatt$assign %in% adrop #columns to drop (always the intercept)
X <- X[, !xdrop, drop=FALSE]
attr(X, "assign") <- Xatt$assign[!xdrop]
if (any(adrop>0)) attr(X, "contrasts") <- Xatt$contrasts[-adrop]
else attr(X, "contrasts") <- Xatt$contrasts
attr(X, "contrasts") <- Xatt$contrasts
offset <- model.offset(mf)
if (is.null(offset) | all(offset==0)) offset <- rep(0., nrow(mf))
assign <- attrassign(X, Terms)
contr.save <- attr(X, "contrasts")
nvars <- ncol(X)
if (missing(init)) init <- rep(0, nvars)
surv.dat <- vector(mode = "list", length = 3)
names(surv.dat) <- c("survival", "X", "Z")
surv.dat[[1]] <- data.frame(Y[,1], Y[,2])
if (is.null(dim(instrument)))
{
instrument <- matrix(instrument, ncol = 1)
}
stopifnot(nrow(X) == nrow(instrument))
colnames(surv.dat[[1]]) <- c("log.t", "delta")
surv.dat[[2]] <- X
surv.dat[[3]] <- instrument
attr(surv.dat, "class") <- "survival.data"
beta <- se.hat <- array(0, dim = c(length(method), nvars))
rownames(beta) <- rownames(se.hat) <- method
colnames(beta) <- colnames(se.hat) <- colnames(X)
fit.objects <- bootstrap.objects <- vector(length = length(method), mode = "list")
names(fit.objects) <- names(bootstrap.objects) <- method
quiet <- TRUE
trace <- FALSE
if (verbose > 1)
{
trace <- TRUE
quiet <- FALSE
}
if (is.null(BB.control))
{
BB.ns <- TRUE
} else
{
BB.ns <- FALSE
}
# solve each estimating equation for beta
for (e in 1:length(method))
{
# return correct estimating equation function
est.eqn <- match.fun(funcs[[match(method[e], types)]])
attr(est.eqn, "name") <- funcs[[match(method[e], types)]]
# and return its smooth counterpart
est.eqn.sm <- match.fun(funcs.sm[[match(method[e], types)]])
attr(est.eqn.sm, "name") <- funcs.sm[[match(method[e], types)]]
if (BB.ns)
{
dfsane.tol <- 0.000001
if (method[e] == "AFT-IPCW")
{
dfsane.tol <- 8
}
BB.control <- list(tol = dfsane.tol,
maxit = 1000,
trace = trace,
M = c(500))
}
ssf <- 1e10
ct <- 0
if (e == 1)
{
init.par <- init
} else
{
init.par <- est$par
#init.par <- NULL
}
if (verbose)
{
cat("Fitting model:", method[e], "\n")
}
# solve for beta using deriv-free spectral method
est <- repFitAFT(tol = dfsane.tol,
data = surv.dat,
est.eqn = est.eqn,
est.eqn.sm = est.eqn.sm,
instrument.names = NULL,
confounded.x.names = confounded.x.names,
maxit = maxit,
dependent.censoring = dependent.censoring,
fit.method = "dfsane",
init.par = init.par,
final.fit = smoothed,
method = c(2),
control = BB.control,
quiet = quiet)
if (bootstrap)
{
if (verbose)
{
cat("Bootstrapping model:", method[e], "\n")
}
# bootstrap estimate
bootstrap.objects[[e]] <- bootstrapVar(est, surv.dat, B = B, method = boot.method,
dependent.censoring = dependent.censoring)
se.hat[e, ] <- bootstrap.objects[[e]]$se.hat
}
beta[e, ] <- est$par
fit.objects[[e]] <- est
}
ret <- list(beta = beta,
se = se.hat,
fit.objects = fit.objects,
bootstrap.objects = bootstrap.objects,
B = B,
Call = Call)
class(ret) <- "aftfits"
ret
}
#' @export
repFitAFT <- function(tol = 5,
maxit = 25,
data,
est.eqn = NULL,
est.eqn.sm = NULL,
instrument.names,
confounded.x.names,
init.par = NULL,
init.method = c("lm", "bisection"),
final.fit = TRUE,
dependent.censoring = FALSE,
fit.method = c("dfsane", "multiStart", "nleqslv", "sane"),
...)
{
# repeatedly calls fitAFT until convergence. decreasing noise
# is added to coefficients at each call to encourage solution
ct <- 0
ssf <- best.ssf <- 1e10
while (ssf > tol & ct <= maxit)
{
ct <- ct + 1
old.ssf <- ssf
#solve for beta using deriv-free spectral method
est <- fitAFT(data = data,
est.eqn = est.eqn,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names,
dependent.censoring = dependent.censoring,
fit.method = fit.method,
init.par = init.par,
...)
ssf <- est$sum.sq.fval
sd <- 2e-3 * min(sqrt(best.ssf), 50)
if (ssf < best.ssf)
{
best.est <- est
best.ssf <- ssf
init.par <- best.est$par + rnorm(length(est$par), sd = sd)
} else
{
init.par <- best.est$par + rnorm(length(est$par), sd = sd)
}
print (sprintf("Current ssf: %g Best ssf: %g, sd: %g", ssf, best.ssf, sd))
}
if (final.fit)
{
best.est <- fitAFT(data = data,
est.eqn = est.eqn.sm,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names,
fit.method = "nleqslv",
init.par = init.par,
global = "dbldog",
method = "Broyden",
dependent.censoring = dependent.censoring,
control = list(ftol = 1e-3,
btol = 1e-4,
trace = 1))
}
best.est$est.eqn <- est.eqn
best.est$est.eqn.smooth <- est.eqn.sm
best.est$final.fit <- final.fit
best.est$call <- match.call()
best.est
}
#' @export
fitAFT <- function(data,
est.eqn = NULL,
instrument.names,
confounded.x.names,
init.par = NULL,
init.method = c("lm", "bisection"),
fit.method = c("dfsane", "multiStart", "nleqslv", "sane"),
dependent.censoring = FALSE,
...)
{
fit.method <- match.arg(fit.method)
if (fit.method == "dfsane" & is.null(est.eqn))
{
stop("Please supply est.eqn if using fit.method: df.sane")
}
GC <- NULL
init.method <- match.arg(init.method)
#instr.loc <- match(instrument.names, colnames(data$X))
conf.x.loc <- match(confounded.x.names, colnames(data$X))
ZXmat <- data$X
if (attr(est.eqn, "name") == "AFT2SLSScorePre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothPre"
| attr(est.eqn, "name") == "AFT2SLSScoreAllPre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothAllPre")
{
#if 2SLS is used, replace Z with Xhat
Z <- as.matrix(data$Z)
if (attr(est.eqn, "name") == "AFT2SLSScoreAllPre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothAllPre")
{
for (v in 1:ncol(data$X))
{
data$X[,v] <- lm(data$X[,v] ~ Z)$fitted.values
}
} else
{
for (v in 1:length(conf.x.loc))
{
dat.lm <- data.frame(response = data$X[,conf.x.loc[v]], predictor = Z)
Xhat <- lm(response ~ predictor, data = dat.lm)$fitted.values
data$X[,conf.x.loc[v]] <- Xhat
}
}
if (attr(est.eqn, "name") == "AFT2SLSScoreSmoothPre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothAllPre")
{
est.eqn <- AFTScoreSmoothPre
attr(est.eqn, "name") <- "AFTScoreSmoothPre"
} else
{
est.eqn <- AFTScorePre
attr(est.eqn, "name") <- "AFTScorePre"
}
} else
{
ZXmat[,conf.x.loc] <- data$Z
}
if (is.null(init.par))
{
#Xhat <- lm(data$X[,conf.x.loc] ~ data$Z)$fitted.values
XXhat <- as.matrix(data$X)
#XXhat[,conf.x.loc] <- Xhat
init.par <- lm(data$survival$log.t ~ XXhat-1)$coefficients
init.par[which(is.na(init.par))] <- 0
names(init.par) <- colnames(data$X)
#num.vars <- length(init.par)
#names(init.par)[conf.x.loc] <- instrument.names
if (init.method == "bisection")
{
interval <- list()
for (i in 1:num.vars)
{
interval[[i]] <- c(init.par[i] - 1, init.par[i] + 1)
}
init.par <- coordinateBisection(est.eqn,
interval = interval,
num.vars = num.vars,
max.iter = 500,
bisection.tol = 0.25,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
print (init.par)
print (sum(est.eqn(init.par, survival=data$survival, X = as.matrix(data$X), ZXmat = as.matrix(ZXmat)) ^ 2))
}
names(init.par) <- colnames(data$X)
}
if (fit.method == "dfsane")
{
#Derivative-Free Spectral Approach for solving nonlinear systems of equations
#from CRAN package 'BB'
if (attr(est.eqn, "name") == "AFTivIPCWScorePre")
{
#generate function G_c() for ICPW
if (dependent.censoring)
{
GC <- genKMCensoringFunc(data$survival, cox = TRUE, X = as.matrix(cbind(data$X, data$Z)))
} else
{
GC <- genKMCensoringFunc(data$survival)
}
df <- BBsolve(par = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
Z = data$Z,
GC = GC,
conf.x.loc = conf.x.loc)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
Z = data$Z,
GC = GC,
conf.x.loc = conf.x.loc)
} else
{
df <- BBsolve(par = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
}
ret <- list(par = df$par,
fval = fval,
iter = df$iter,
nobs = nrow(data$X),
nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
} else if (fit.method == "nleqslv" | fit.method == "sane")
{
if (!is.null(est.eqn))
{
if (attr(est.eqn, "name") != "AFTScoreSmoothPre" & attr(est.eqn, "name") != "AFTivScoreSmoothPre")
{
warning(paste(paste("Arg: est.eqn =", attr(est.eqn, "name")),
"not be used. Used AFTivScoreSmoothPre instead"))
est.eqn <- AFTivScoreSmoothPre
}
} else
{
warning("est.eqn not supplied. Used AFTivScoreSmoothPre")
est.eqn <- AFTivScoreSmoothPre
}
if (fit.method == "nleqslv")
{
df <- nleqslv(x = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
tau = 0.01)
ret <- list(par = df$x,
fval = df$fvec,
iter = df$iter,
nobs = nrow(data$X),
nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
#print(df)
} else if (fit.method == "sane")
{
df <- sane(par = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
tau = 0.01)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
tau = 0.01)
ret <- list(par = df$par,
fval = fval,
iter = df$iter,
nobs = nrow(data$X),
nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
}
} else if (fit.method == "multiStart")
{
#Derivative-Free Spectral Approach for solving nonlinear systems of equations
#from CRAN package 'BB'
n.starts <- 10
p <- length(init.par)
p0 <- matrix(rnorm(n.starts * p), n.starts, p)
p0 <- rbind(p0, init.par)
if (attr(est.eqn, "name") == "AFTivIPCWScorePre")
{
#generate function G_c() for ICPW
if (dependent.censoring)
{
GC <- genKMCensoringFunc(data$survival, cox = TRUE, X = as.matrix(cbind(data$X, data$Z)))
} else
{
GC <- genKMCensoringFunc(data$survival)
}
df <- multiStart(par = p0,
fn = est.eqn,
...,
action = "solve",
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
Z = as.matrix(data$Z),
GC = GC)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
GC = GC)
} else
{
df <- multiStart(par = p0,
fn = est.eqn,
...,
action = "solve",
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
#take only the best ones and fit again
good.init.idx <- which(df$fvalue < 500)
p0 <- df$par[good.init.idx,] + matrix(rnorm(length(good.init.idx) * p, sd = 5e-6), length(good.init.idx), p)
df <- multiStart(par = p0,
fn = est.eqn,
...,
action = "solve",
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
good.init.idx <- which(df$fvalue < 100)
df$par <- df$par[good.init.idx,]
df$fvalue <- df$fvalue[good.init.idx]
colnames(p0) <- colnames(data$X)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
}
ret <- list(par = df$par, fval = fval, iter = df$iter,
nobs = nrow(data$X), nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
}
ret$sum.sq.fval <- sum(ret$fval^2)
if (init.method == "multiStart")
{
ret <- df
}
ret$call <- match.call()
class(ret) <- "aft.fit"
ret
}
#' @export
bootstrap.model <- function(model,
formula,
data,
instrument,
confounded.x.names,
method = c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW"),
smoothed = FALSE,
weights = NULL,
bootstrap = FALSE,
boot.method = c("ls", "sv", "full.bootstrap"),
B = 100L,
na.action,
init = NULL,
return.data = FALSE,
tol = 1e-5,
maxit = 100,
verbose = 0,
BB.control = NULL,
...)
{
method <- match.arg(method, several.ok = TRUE)
boot.method <- match.arg(boot.method)
Call <- match.call()
# create a call to model.frame() that contains the formula (required)
# and any other of the relevant optional arguments
# then evaluate it in the proper frame
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
names(Call), nomatch = 0)
if (indx[1] == 0) stop("A formula argument is required")
temp <- Call[c(1,indx)] # only keep the arguments we wanted
temp[[1]] <- as.name('model.frame') # change the function called
temp$formula <- if(missing(data)) terms(formula)
else terms(formula, data=data)
mf <- eval(temp, parent.frame())
if (nrow(mf) == 0) stop("No (non-missing) observations")
Terms <- terms(mf)
Y <- model.extract(mf, "response")
if (!inherits(Y, "Surv")) stop("Response must be a survival object")
type <- attr(Y, "type")
if (type!='right' && type!='counting')
{
stop(paste("AFT model doesn't support \"", type,
"\" survival data", sep=''))
}
weights <- model.weights(mf)
nobs <- nrow(Y) #remember this before any time transforms
types <- c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW")
funcs <- c("AFTScorePre", "AFTivScorePre", "AFT2SLSScorePre", "AFTivIPCWScorePre")
funcs.sm <- c("AFTScoreSmoothPre", "AFTivScoreSmoothPre", "AFT2SLSScoreSmoothPre", "AFTivIPCWScoreSmooth")
contrast.arg <- NULL #due to shared code with model.matrix.coxph
X <- model.matrix(Terms, mf, contrasts=contrast.arg)
#print(head(X))
# drop the intercept after the fact, and also drop strata if necessary
Xatt <- attributes(X)
adrop <- 0
xdrop <- Xatt$assign %in% adrop
#xdrop <- Xatt$assign %in% adrop #columns to drop (always the intercept)
X <- X[, !xdrop, drop=FALSE]
attr(X, "assign") <- Xatt$assign[!xdrop]
if (any(adrop>0)) attr(X, "contrasts") <- Xatt$contrasts[-adrop]
else attr(X, "contrasts") <- Xatt$contrasts
attr(X, "contrasts") <- Xatt$contrasts
offset <- model.offset(mf)
if (is.null(offset) | all(offset==0)) offset <- rep(0., nrow(mf))
assign <- attrassign(X, Terms)
contr.save <- attr(X, "contrasts")
nvars <- ncol(X)
if (missing(init)) init <- rep(0, nvars)
surv.dat <- vector(mode = "list", length = 3)
names(surv.dat) <- c("survival", "X", "Z")
surv.dat[[1]] <- data.frame(Y[,1], Y[,2])
if (is.null(dim(instrument)))
{
instrument <- matrix(instrument, ncol = 1)
}
stopifnot(nrow(X) == nrow(instrument))
colnames(surv.dat[[1]]) <- c("log.t", "delta")
surv.dat[[2]] <- X
surv.dat[[3]] <- instrument
attr(surv.dat, "class") <- "survival.data"
beta <- se.hat <- array(0, dim = c(length(method), nvars))
rownames(beta) <- rownames(se.hat) <- method
colnames(beta) <- colnames(se.hat) <- colnames(X)
fit.objects <- bootstrap.objects <- vector(length = length(method), mode = "list")
names(fit.objects) <- names(bootstrap.objects) <- method
quiet <- TRUE
trace <- FALSE
if (verbose > 1)
{
trace <- TRUE
quiet <- FALSE
}
if (is.null(BB.control))
{
BB.ns <- TRUE
} else
{
BB.ns <- FALSE
}
#solve each estimating equation for beta
for (e in 1:length(method))
{
# return correct estimating equation function
est.eqn <- match.fun(funcs[[match(method[e], types)]])
attr(est.eqn, "name") <- funcs[[match(method[e], types)]]
# and return the smooth counterpart
est.eqn.sm <- match.fun(funcs.sm[[match(method[e], types)]])
attr(est.eqn.sm, "name") <- funcs.sm[[match(method[e], types)]]
ct <- 0
if (bootstrap)
{
if (verbose)
{
cat("Bootstrapping model", e, "\n")
}
# bootstrap estimate
bootstrap.objects[[e]] <- bootstrapVar(model$fit.objects[[e]], surv.dat, B = B, method = boot.method)
se.hat[e, ] <- bootstrap.objects[[e]]$se.hat
}
#beta[e, ] <- est$par
#fit.objects[[e]] <- est
}
ret <- model
ret$bootstrap.objects <- bootstrap.objects
ret$se <- se.hat
class(ret) <- "aftiv.boot.fit"
ret
}
jaredhuling/aftiv documentation built on May 20, 2019, 9:55 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions aftfit repFitAFT fitAFT bootstrap.model
Documented in aftfit
#' Instrumental variable estimation via semiparametric rank-based AFT models
#'
#' @description Fitting function for instrumental variable estimation under the semiparametric AFT model
#' @param formula a formula object, with the response on the left of a tilde operator, and the covariates.
#' The response must be a survival object as returned by the \code{Surv()} function.
#' @param data a \code{data.frame} which contains the variables named in the \code{formula}.
#' @param instrument a vector of length equal to the number of rows in \code{data} containing the instrumental variable
#' @param confounded.x.names name of the exposure variable (endogenous variable) of interest
#' @param method one of:
#' \itemize{
#' \item{\code{"AFT-IPCW"}}{ The proposed IV method that uses inverse probability of censoring weighting }
#' \item{\code{"AFT"}} { Plain rank-based AFT method that does not account for unmeasured confounding}
#' \item{\code{"AFT-IV"}} { biased AFT method that acccounts for unmeasured confounding, but not the non-ignorable censoring
#' (and is hence biased) }
#' \item{\code{"AFT-2SLS"}} { Two-stage AFT method which accounts for unmeasured confounding but requires correct specification
#' of relationship between IV and exposure }
#' }
#' @param smoothed should a smoothed version of the estimating equation be used? Experimental.
#' @param weights observation weights. Not used currently!
#' @param bootstrap logical variable indicating whether or not to run bootstrap
#' @param boot.method bootstrap method to use. \code{"ls"} uses the LS method of Zeng and Lin (2008),
#' \code{"sv"} uses the SV method of Zeng and Lin (2008). \code{"full.bootstrap"} uses a full boostrap where resampled estimating equations
#' are solved for each iteration (this method is very slow!)
#' @param B number of bootstrap iterations
#' @param dependent.censoring for the \code{"AFT-IPCW"} method, should the censoring model allow for dependence on covariates? If \code{TRUE}
#' a Cox model will be fit for the censoring distribution
#' @param na.action how should missing data be treated? See help file for \code{\link[stats]{lm}} for more details
#' @param return.data should the data be returned?
#' @param init vector equal to the length of the number of parameters in the model used to initialize estimation
#' @param tol positive tolerance threshold. Smaller values indicate more precise solutions are required
#' @param maxit maximum number of restarts for \code{BB} algorithm
#' @param verbose should messages be printed? logical variable
#' @param BB.control control list for \code{BBsolve()} function. See \code{\link[BB]{BBsolve}} for more information
#' @param ... not used
#' @import ggplot2
#' @import stats
#' @export
#'
#'
#' @examples
#' library(aftiv)
#'
#' set.seed(1)
#' true.beta <- c(1,-0.5,-0.5,0.25,0,0,-0.75,0.75)
#' dat <- simIVMultivarSurvivalData(500,1,1,-1,1,true.beta,num.confounded = 1,
#' confounding.function = "exp")
#'
#' df <- data.frame(dat$survival[c("delta", "log.t")], dat$X)
#' Z <- dat$Z
#'
#' system.time(aftf <- aftfit(Surv(log.t, delta) ~ ., data = df,
#' instrument = Z,
#' confounded.x.names = "X1",
#' method = c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW"),
#' boot.method = "ls", verbose = 1,
#' B = 200L, smoothed = FALSE,
#' bootstrap = TRUE))
#'
#' summary(aftf)
aftfit <- function(formula,
data,
instrument,
confounded.x.names,
method = c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW"),
smoothed = FALSE,
weights = NULL,
bootstrap = FALSE,
boot.method = c("ls", "sv", "full.bootstrap"),
B = 1000L,
dependent.censoring = FALSE,
na.action,
init = NULL,
return.data = FALSE,
tol = 1e-5,
maxit = 10L,
verbose = 0,
BB.control = NULL,
...) {
method <- match.arg(method, several.ok = TRUE)
boot.method <- match.arg(boot.method)
Call <- match.call()
# create a call to model.frame() that contains the formula (required)
# and any other of the relevant optional arguments
# then evaluate it in the proper frame
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
names(Call), nomatch=0)
if (indx[1] == 0) stop("A formula argument is required")
temp <- Call[c(1,indx)] # only keep the arguments we wanted
temp[[1]] <- as.name('model.frame') # change the function called
temp$formula <- if(missing(data)) terms(formula)
else terms(formula, data=data)
mf <- eval(temp, parent.frame())
if (nrow(mf) == 0) stop("No (non-missing) observations")
Terms <- terms(mf)
Y <- model.extract(mf, "response")
if (!inherits(Y, "Surv")) stop("Response must be a survival object")
type <- attr(Y, "type")
if (type!='right' && type!='counting')
stop(paste("AFT model doesn't support \"", type,
"\" survival data", sep=''))
weights <- model.weights(mf)
nobs <- nrow(Y) #remember this before any time transforms
types <- c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW")
funcs <- c("AFTScorePre", "AFTivScorePre", "AFT2SLSScorePre", "AFTivIPCWScorePre")
funcs.sm <- c("AFTScoreSmoothPre", "AFTivScoreSmoothPre", "AFT2SLSScoreSmoothPre", "AFTivIPCWScoreSmooth")
contrast.arg <- NULL #due to shared code with model.matrix.coxph
X <- model.matrix(Terms, mf, contrasts=contrast.arg)
#print(head(X))
# drop the intercept after the fact, and also drop strata if necessary
Xatt <- attributes(X)
adrop <- 0
xdrop <- Xatt$assign %in% adrop
#xdrop <- Xatt$assign %in% adrop #columns to drop (always the intercept)
X <- X[, !xdrop, drop=FALSE]
attr(X, "assign") <- Xatt$assign[!xdrop]
if (any(adrop>0)) attr(X, "contrasts") <- Xatt$contrasts[-adrop]
else attr(X, "contrasts") <- Xatt$contrasts
attr(X, "contrasts") <- Xatt$contrasts
offset <- model.offset(mf)
if (is.null(offset) | all(offset==0)) offset <- rep(0., nrow(mf))
assign <- attrassign(X, Terms)
contr.save <- attr(X, "contrasts")
nvars <- ncol(X)
if (missing(init)) init <- rep(0, nvars)
surv.dat <- vector(mode = "list", length = 3)
names(surv.dat) <- c("survival", "X", "Z")
surv.dat[[1]] <- data.frame(Y[,1], Y[,2])
if (is.null(dim(instrument)))
{
instrument <- matrix(instrument, ncol = 1)
}
stopifnot(nrow(X) == nrow(instrument))
colnames(surv.dat[[1]]) <- c("log.t", "delta")
surv.dat[[2]] <- X
surv.dat[[3]] <- instrument
attr(surv.dat, "class") <- "survival.data"
beta <- se.hat <- array(0, dim = c(length(method), nvars))
rownames(beta) <- rownames(se.hat) <- method
colnames(beta) <- colnames(se.hat) <- colnames(X)
fit.objects <- bootstrap.objects <- vector(length = length(method), mode = "list")
names(fit.objects) <- names(bootstrap.objects) <- method
quiet <- TRUE
trace <- FALSE
if (verbose > 1)
{
trace <- TRUE
quiet <- FALSE
}
if (is.null(BB.control))
{
BB.ns <- TRUE
} else
{
BB.ns <- FALSE
}
# solve each estimating equation for beta
for (e in 1:length(method))
{
# return correct estimating equation function
est.eqn <- match.fun(funcs[[match(method[e], types)]])
attr(est.eqn, "name") <- funcs[[match(method[e], types)]]
# and return its smooth counterpart
est.eqn.sm <- match.fun(funcs.sm[[match(method[e], types)]])
attr(est.eqn.sm, "name") <- funcs.sm[[match(method[e], types)]]
if (BB.ns)
{
dfsane.tol <- 0.000001
if (method[e] == "AFT-IPCW")
{
dfsane.tol <- 8
}
BB.control <- list(tol = dfsane.tol,
maxit = 1000,
trace = trace,
M = c(500))
}
ssf <- 1e10
ct <- 0
if (e == 1)
{
init.par <- init
} else
{
init.par <- est$par
#init.par <- NULL
}
if (verbose)
{
cat("Fitting model:", method[e], "\n")
}
# solve for beta using deriv-free spectral method
est <- repFitAFT(tol = dfsane.tol,
data = surv.dat,
est.eqn = est.eqn,
est.eqn.sm = est.eqn.sm,
instrument.names = NULL,
confounded.x.names = confounded.x.names,
maxit = maxit,
dependent.censoring = dependent.censoring,
fit.method = "dfsane",
init.par = init.par,
final.fit = smoothed,
method = c(2),
control = BB.control,
quiet = quiet)
if (bootstrap)
{
if (verbose)
{
cat("Bootstrapping model:", method[e], "\n")
}
# bootstrap estimate
bootstrap.objects[[e]] <- bootstrapVar(est, surv.dat, B = B, method = boot.method,
dependent.censoring = dependent.censoring)
se.hat[e, ] <- bootstrap.objects[[e]]$se.hat
}
beta[e, ] <- est$par
fit.objects[[e]] <- est
}
ret <- list(beta = beta,
se = se.hat,
fit.objects = fit.objects,
bootstrap.objects = bootstrap.objects,
B = B,
Call = Call)
class(ret) <- "aftfits"
ret
}
#' @export
repFitAFT <- function(tol = 5,
maxit = 25,
data,
est.eqn = NULL,
est.eqn.sm = NULL,
instrument.names,
confounded.x.names,
init.par = NULL,
init.method = c("lm", "bisection"),
final.fit = TRUE,
dependent.censoring = FALSE,
fit.method = c("dfsane", "multiStart", "nleqslv", "sane"),
...)
{
# repeatedly calls fitAFT until convergence. decreasing noise
# is added to coefficients at each call to encourage solution
ct <- 0
ssf <- best.ssf <- 1e10
while (ssf > tol & ct <= maxit)
{
ct <- ct + 1
old.ssf <- ssf
#solve for beta using deriv-free spectral method
est <- fitAFT(data = data,
est.eqn = est.eqn,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names,
dependent.censoring = dependent.censoring,
fit.method = fit.method,
init.par = init.par,
...)
ssf <- est$sum.sq.fval
sd <- 2e-3 * min(sqrt(best.ssf), 50)
if (ssf < best.ssf)
{
best.est <- est
best.ssf <- ssf
init.par <- best.est$par + rnorm(length(est$par), sd = sd)
} else
{
init.par <- best.est$par + rnorm(length(est$par), sd = sd)
}
print (sprintf("Current ssf: %g Best ssf: %g, sd: %g", ssf, best.ssf, sd))
}
if (final.fit)
{
best.est <- fitAFT(data = data,
est.eqn = est.eqn.sm,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names,
fit.method = "nleqslv",
init.par = init.par,
global = "dbldog",
method = "Broyden",
dependent.censoring = dependent.censoring,
control = list(ftol = 1e-3,
btol = 1e-4,
trace = 1))
}
best.est$est.eqn <- est.eqn
best.est$est.eqn.smooth <- est.eqn.sm
best.est$final.fit <- final.fit
best.est$call <- match.call()
best.est
}
#' @export
fitAFT <- function(data,
est.eqn = NULL,
instrument.names,
confounded.x.names,
init.par = NULL,
init.method = c("lm", "bisection"),
fit.method = c("dfsane", "multiStart", "nleqslv", "sane"),
dependent.censoring = FALSE,
...)
{
fit.method <- match.arg(fit.method)
if (fit.method == "dfsane" & is.null(est.eqn))
{
stop("Please supply est.eqn if using fit.method: df.sane")
}
GC <- NULL
init.method <- match.arg(init.method)
#instr.loc <- match(instrument.names, colnames(data$X))
conf.x.loc <- match(confounded.x.names, colnames(data$X))
ZXmat <- data$X
if (attr(est.eqn, "name") == "AFT2SLSScorePre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothPre"
| attr(est.eqn, "name") == "AFT2SLSScoreAllPre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothAllPre")
{
#if 2SLS is used, replace Z with Xhat
Z <- as.matrix(data$Z)
if (attr(est.eqn, "name") == "AFT2SLSScoreAllPre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothAllPre")
{
for (v in 1:ncol(data$X))
{
data$X[,v] <- lm(data$X[,v] ~ Z)$fitted.values
}
} else
{
for (v in 1:length(conf.x.loc))
{
dat.lm <- data.frame(response = data$X[,conf.x.loc[v]], predictor = Z)
Xhat <- lm(response ~ predictor, data = dat.lm)$fitted.values
data$X[,conf.x.loc[v]] <- Xhat
}
}
if (attr(est.eqn, "name") == "AFT2SLSScoreSmoothPre" | attr(est.eqn, "name") == "AFT2SLSScoreSmoothAllPre")
{
est.eqn <- AFTScoreSmoothPre
attr(est.eqn, "name") <- "AFTScoreSmoothPre"
} else
{
est.eqn <- AFTScorePre
attr(est.eqn, "name") <- "AFTScorePre"
}
} else
{
ZXmat[,conf.x.loc] <- data$Z
}
if (is.null(init.par))
{
#Xhat <- lm(data$X[,conf.x.loc] ~ data$Z)$fitted.values
XXhat <- as.matrix(data$X)
#XXhat[,conf.x.loc] <- Xhat
init.par <- lm(data$survival$log.t ~ XXhat-1)$coefficients
init.par[which(is.na(init.par))] <- 0
names(init.par) <- colnames(data$X)
#num.vars <- length(init.par)
#names(init.par)[conf.x.loc] <- instrument.names
if (init.method == "bisection")
{
interval <- list()
for (i in 1:num.vars)
{
interval[[i]] <- c(init.par[i] - 1, init.par[i] + 1)
}
init.par <- coordinateBisection(est.eqn,
interval = interval,
num.vars = num.vars,
max.iter = 500,
bisection.tol = 0.25,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
print (init.par)
print (sum(est.eqn(init.par, survival=data$survival, X = as.matrix(data$X), ZXmat = as.matrix(ZXmat)) ^ 2))
}
names(init.par) <- colnames(data$X)
}
if (fit.method == "dfsane")
{
#Derivative-Free Spectral Approach for solving nonlinear systems of equations
#from CRAN package 'BB'
if (attr(est.eqn, "name") == "AFTivIPCWScorePre")
{
#generate function G_c() for ICPW
if (dependent.censoring)
{
GC <- genKMCensoringFunc(data$survival, cox = TRUE, X = as.matrix(cbind(data$X, data$Z)))
} else
{
GC <- genKMCensoringFunc(data$survival)
}
df <- BBsolve(par = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
Z = data$Z,
GC = GC,
conf.x.loc = conf.x.loc)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
Z = data$Z,
GC = GC,
conf.x.loc = conf.x.loc)
} else
{
df <- BBsolve(par = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
}
ret <- list(par = df$par,
fval = fval,
iter = df$iter,
nobs = nrow(data$X),
nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
} else if (fit.method == "nleqslv" | fit.method == "sane")
{
if (!is.null(est.eqn))
{
if (attr(est.eqn, "name") != "AFTScoreSmoothPre" & attr(est.eqn, "name") != "AFTivScoreSmoothPre")
{
warning(paste(paste("Arg: est.eqn =", attr(est.eqn, "name")),
"not be used. Used AFTivScoreSmoothPre instead"))
est.eqn <- AFTivScoreSmoothPre
}
} else
{
warning("est.eqn not supplied. Used AFTivScoreSmoothPre")
est.eqn <- AFTivScoreSmoothPre
}
if (fit.method == "nleqslv")
{
df <- nleqslv(x = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
tau = 0.01)
ret <- list(par = df$x,
fval = df$fvec,
iter = df$iter,
nobs = nrow(data$X),
nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
#print(df)
} else if (fit.method == "sane")
{
df <- sane(par = init.par,
fn = est.eqn,
...,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
tau = 0.01)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
tau = 0.01)
ret <- list(par = df$par,
fval = fval,
iter = df$iter,
nobs = nrow(data$X),
nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
}
} else if (fit.method == "multiStart")
{
#Derivative-Free Spectral Approach for solving nonlinear systems of equations
#from CRAN package 'BB'
n.starts <- 10
p <- length(init.par)
p0 <- matrix(rnorm(n.starts * p), n.starts, p)
p0 <- rbind(p0, init.par)
if (attr(est.eqn, "name") == "AFTivIPCWScorePre")
{
#generate function G_c() for ICPW
if (dependent.censoring)
{
GC <- genKMCensoringFunc(data$survival, cox = TRUE, X = as.matrix(cbind(data$X, data$Z)))
} else
{
GC <- genKMCensoringFunc(data$survival)
}
df <- multiStart(par = p0,
fn = est.eqn,
...,
action = "solve",
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
Z = as.matrix(data$Z),
GC = GC)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat),
GC = GC)
} else
{
df <- multiStart(par = p0,
fn = est.eqn,
...,
action = "solve",
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
#take only the best ones and fit again
good.init.idx <- which(df$fvalue < 500)
p0 <- df$par[good.init.idx,] + matrix(rnorm(length(good.init.idx) * p, sd = 5e-6), length(good.init.idx), p)
df <- multiStart(par = p0,
fn = est.eqn,
...,
action = "solve",
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
good.init.idx <- which(df$fvalue < 100)
df$par <- df$par[good.init.idx,]
df$fvalue <- df$fvalue[good.init.idx]
colnames(p0) <- colnames(data$X)
fval <- est.eqn(beta = df$par,
survival = data$survival,
X = as.matrix(data$X),
ZXmat = as.matrix(ZXmat))
}
ret <- list(par = df$par, fval = fval, iter = df$iter,
nobs = nrow(data$X), nvars = length(df$par),
GC = GC,
instrument.names = instrument.names,
confounded.x.names = confounded.x.names)
}
ret$sum.sq.fval <- sum(ret$fval^2)
if (init.method == "multiStart")
{
ret <- df
}
ret$call <- match.call()
class(ret) <- "aft.fit"
ret
}
#' @export
bootstrap.model <- function(model,
formula,
data,
instrument,
confounded.x.names,
method = c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW"),
smoothed = FALSE,
weights = NULL,
bootstrap = FALSE,
boot.method = c("ls", "sv", "full.bootstrap"),
B = 100L,
na.action,
init = NULL,
return.data = FALSE,
tol = 1e-5,
maxit = 100,
verbose = 0,
BB.control = NULL,
...)
{
method <- match.arg(method, several.ok = TRUE)
boot.method <- match.arg(boot.method)
Call <- match.call()
# create a call to model.frame() that contains the formula (required)
# and any other of the relevant optional arguments
# then evaluate it in the proper frame
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
names(Call), nomatch = 0)
if (indx[1] == 0) stop("A formula argument is required")
temp <- Call[c(1,indx)] # only keep the arguments we wanted
temp[[1]] <- as.name('model.frame') # change the function called
temp$formula <- if(missing(data)) terms(formula)
else terms(formula, data=data)
mf <- eval(temp, parent.frame())
if (nrow(mf) == 0) stop("No (non-missing) observations")
Terms <- terms(mf)
Y <- model.extract(mf, "response")
if (!inherits(Y, "Surv")) stop("Response must be a survival object")
type <- attr(Y, "type")
if (type!='right' && type!='counting')
{
stop(paste("AFT model doesn't support \"", type,
"\" survival data", sep=''))
}
weights <- model.weights(mf)
nobs <- nrow(Y) #remember this before any time transforms
types <- c("AFT", "AFT-IV", "AFT-2SLS", "AFT-IPCW")
funcs <- c("AFTScorePre", "AFTivScorePre", "AFT2SLSScorePre", "AFTivIPCWScorePre")
funcs.sm <- c("AFTScoreSmoothPre", "AFTivScoreSmoothPre", "AFT2SLSScoreSmoothPre", "AFTivIPCWScoreSmooth")
contrast.arg <- NULL #due to shared code with model.matrix.coxph
X <- model.matrix(Terms, mf, contrasts=contrast.arg)
#print(head(X))
# drop the intercept after the fact, and also drop strata if necessary
Xatt <- attributes(X)
adrop <- 0
xdrop <- Xatt$assign %in% adrop
#xdrop <- Xatt$assign %in% adrop #columns to drop (always the intercept)
X <- X[, !xdrop, drop=FALSE]
attr(X, "assign") <- Xatt$assign[!xdrop]
if (any(adrop>0)) attr(X, "contrasts") <- Xatt$contrasts[-adrop]
else attr(X, "contrasts") <- Xatt$contrasts
attr(X, "contrasts") <- Xatt$contrasts
offset <- model.offset(mf)
if (is.null(offset) | all(offset==0)) offset <- rep(0., nrow(mf))
assign <- attrassign(X, Terms)
contr.save <- attr(X, "contrasts")
nvars <- ncol(X)
if (missing(init)) init <- rep(0, nvars)
surv.dat <- vector(mode = "list", length = 3)
names(surv.dat) <- c("survival", "X", "Z")
surv.dat[[1]] <- data.frame(Y[,1], Y[,2])
if (is.null(dim(instrument)))
{
instrument <- matrix(instrument, ncol = 1)
}
stopifnot(nrow(X) == nrow(instrument))
colnames(surv.dat[[1]]) <- c("log.t", "delta")
surv.dat[[2]] <- X
surv.dat[[3]] <- instrument
attr(surv.dat, "class") <- "survival.data"
beta <- se.hat <- array(0, dim = c(length(method), nvars))
rownames(beta) <- rownames(se.hat) <- method
colnames(beta) <- colnames(se.hat) <- colnames(X)
fit.objects <- bootstrap.objects <- vector(length = length(method), mode = "list")
names(fit.objects) <- names(bootstrap.objects) <- method
quiet <- TRUE
trace <- FALSE
if (verbose > 1)
{
trace <- TRUE
quiet <- FALSE
}
if (is.null(BB.control))
{
BB.ns <- TRUE
} else
{
BB.ns <- FALSE
}
#solve each estimating equation for beta
for (e in 1:length(method))
{
# return correct estimating equation function
est.eqn <- match.fun(funcs[[match(method[e], types)]])
attr(est.eqn, "name") <- funcs[[match(method[e], types)]]
# and return the smooth counterpart
est.eqn.sm <- match.fun(funcs.sm[[match(method[e], types)]])
attr(est.eqn.sm, "name") <- funcs.sm[[match(method[e], types)]]
ct <- 0
if (bootstrap)
{
if (verbose)
{
cat("Bootstrapping model", e, "\n")
}
# bootstrap estimate
bootstrap.objects[[e]] <- bootstrapVar(model$fit.objects[[e]], surv.dat, B = B, method = boot.method)
se.hat[e, ] <- bootstrap.objects[[e]]$se.hat
}
#beta[e, ] <- est$par
#fit.objects[[e]] <- est
}
ret <- model
ret$bootstrap.objects <- bootstrap.objects
ret$se <- se.hat
class(ret) <- "aftiv.boot.fit"
ret
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.