Nothing
R/deriv_functions.R
In CIEE: Estimating and Testing Direct Effects in Directed Acyclic Graphs using Estimating Equations
Defines functions
deriv_obj
scores
hessian
Documented in
deriv_obj
hessian
scores
#' Score and hessian matrix based on the estimating functions.
#'
#' Functions to compute the score and hessian matrices of the parameters
#' based on the estimating functions, under the GLM and AFT setting for
#' the analysis of a normally-distributed or censored time-to-event
#' primary outcome. The score and hessian matrices are further used in
#' the functions \code{\link{sandwich_se}}, \code{\link{ciee}} and
#' \code{\link{ciee_loop}} to obtain robust sandwich error estimates of the
#' parameter estimates of
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1^2, \alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1^2, \alpha4, \alphaXY, \sigma2^2}
#' under the GLM setting and
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1, \alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1, \alpha4, \alphaXY, \sigma2^2}
#' under the AFT setting.
#'
#' For the computation of the score and hessian matrices, first, the help function
#' \code{\link{deriv_obj}} is used. In a first step, the expression of all first
#' and second derivatives of the parameters is computed using the expressions of
#' \code{logL1} and \code{logL2} from the \code{\link{est_funct_expr}} as input.
#' Then, the numerical values of all first and second derivatives are obtained
#' for the observed data \code{Y}, \code{X}, \code{K}, \code{L} (and \code{C} under
#' the AFT setting) and point estimates (\code{estimates}) of the parameters,
#' for all observed individuals.
#'
#' Second, the functions \code{\link{scores}} and \code{\link{hessian}} are used
#' to extract the relevant score and hessian matrices with respect to \code{logL1}
#' and \code{logL2} from the output of \code{\link{deriv_obj}} and piece them together.
#' For further details, see the vignette.
#'
#' @param setting String with value \code{"GLM"} or \code{"AFT"} indicating
#' whether the matrices are computed under the GLM or AFT setting.
#' @param logL1 Expression of the function \code{logL1} generated by
#' the \code{\link{est_funct_expr}} function.
#' @param logL2 Expression of the function \code{logL2} generated by
#' the \code{\link{est_funct_expr}} function.
#' @param Y Numeric input vector for the primary outcome.
#' @param X Numeric input vector for the exposure variable.
#' @param K Numeric input vector for the intermediate outcome.
#' @param L Numeric input vector for the observed confounding factor.
#' @param C Numeric input vector for the censoring indicator under the AFT setting
#' (must be coded 0 = censored, 1 = uncensored).
#' @param estimates Numeric input vector with point estimates of the parameters
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1^2,}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1^2,}
#' \eqn{\alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha4, \alphaXY, \sigma2^2}
#' under the GLM setting and of
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1, \alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1, \alpha4, \alphaXY, \sigma2^2}
#' under the AFT setting. Under the AFT setting,
#' \code{estimates} must also contain the mean of the estimated
#' true survival times \code{"y_adj_bar"}.
#' @param derivobj Output of the \code{\link{deriv_obj}} function used as input in
#' the \code{\link{scores}} and \code{\link{hessian}} functions.
#'
#' @return The \code{\link{deriv_obj}} function returns a list with
#' objects \code{logL1_deriv}, \code{logL2_deriv} which
#' contain the score and hessian matrices based on \code{logL1},
#' \code{logL2}, respectively.
#' @return The \code{\link{scores}} function returns the \eqn{(n \times 8)}{(n x 8)}
#' score matrix.
#' @return The \code{\link{hessian}} function returns the \eqn{(n \times 8 \times 8)}{(n x 8 x 8)}
#' hessian matrix.
#'
#' @examples
#'
#' # Generate data including Y, K, L, X under the GLM setting
#' dat <- generate_data(setting = "GLM")
#'
#' # Obtain estimating functions' expressions
#' estfunct <- est_funct_expr(setting = "GLM")
#'
#' # Obtain point estimates of the parameters
#' estimates <- get_estimates(setting = "GLM", Y = dat$Y, X = dat$X,
#' K = dat$K, L = dat$L)
#'
#' # Obtain matrices with all first and second derivatives
#' derivobj <- deriv_obj(setting = "GLM", logL1 = estfunct$logL1,
#' logL2 = estfunct$logL2, Y = dat$Y, X = dat$X,
#' K = dat$K, L = dat$L, estimates = estimates)
#' names(derivobj)
#' head(derivobj$logL1_deriv$gradient)
#'
#' # Obtain score and hessian matrices
#' scores(derivobj)
#' hessian(derivobj)
#'
#' @name score_and_hessian_matrix_functions
NULL
#' @rdname score_and_hessian_matrix_functions
#' @export
deriv_obj <- function(setting = "GLM", logL1 = NULL, logL2 = NULL, Y = NULL,
X = NULL, K = NULL, L = NULL, C = NULL,
estimates = NULL) {
if (is.null(setting) | is.null(logL1) | is.null(logL2) |
is.null(estimates)) {
stop("One or more arguments of the function are missing.")
}
if (is.null(Y) | is.null(X) | is.null(K) | is.null(L)) {
stop("Data of one or more variables are not supplied.")
}
if (setting == "AFT" & is.null(C)) {
stop("C has to be supplied for the AFT setting.")
}
n <- length(Y)
if (setting == "GLM") {
data_help <- data.frame(Y = Y, X = X, K = K, L = L)
data_help <- data_help[stats::complete.cases(data_help), ]
U12345_i <- stats::deriv(expr = logL1, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1sq", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "k_i", "x_i", "l_i", "alpha0",
"alpha1", "alpha2", "alpha3", "sigma1sq"), func = T,
hessian = T)
U678_i <- stats::deriv(expr = logL2, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1sq", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "k_i", "x_i", "y_bar", "k_bar",
"alpha1", "alpha4", "alphaXY", "sigma2sq"), func = T,
hessian = T)
logL1_deriv <- attributes(U12345_i(y_i = data_help$Y, k_i = data_help$K,
x_i = data_help$X, l_i = data_help$L,
alpha0 = estimates[names(estimates) == "alpha_0"],
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha2 = estimates[names(estimates) == "alpha_2"],
alpha3 = estimates[names(estimates) == "alpha_3"],
sigma1sq = estimates[names(estimates) == "sigma_1_sq"]))
logL2_deriv <- attributes(U678_i(y_i = data_help$Y, k_i = data_help$K,
x_i = data_help$X, y_bar = mean(data_help$Y),
k_bar = mean(data_help$K),
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha4 = estimates[names(estimates) == "alpha_4"],
alphaXY = estimates[names(estimates) == "alpha_XY"],
sigma2sq = estimates[names(estimates) == "sigma_2_sq"]))
deriv_obj <- list(logL1_deriv = logL1_deriv, logL2_deriv = logL2_deriv)
}
if (setting == "AFT") {
data_help <- data.frame(Y = Y, X = X, K = K, L = L, C = C)
data_help <- data_help[stats::complete.cases(data_help), ]
U12345_i <- stats::deriv(expr = logL1, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "c_i", "k_i", "x_i", "l_i",
"alpha0", "alpha1", "alpha2", "alpha3", "sigma1"),
func = T, hessian = T)
U678_i <- stats::deriv(expr = logL2, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "c_i", "k_i", "x_i", "l_i",
"y_adj_bar", "k_bar", "alpha0", "alpha1", "alpha2",
"alpha3", "sigma1", "alpha4", "alphaXY", "sigma2sq"),
func = T, hessian = T)
logL1_deriv <- attributes(U12345_i(y_i = data_help$Y, c_i = data_help$C,
k_i = data_help$K, x_i = data_help$X, l_i = data_help$L,
alpha0 = estimates[names(estimates) == "alpha_0"],
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha2 = estimates[names(estimates) == "alpha_2"],
alpha3 = estimates[names(estimates) == "alpha_3"],
sigma1 = estimates[names(estimates) == "sigma_1"]))
logL2_deriv <- attributes(U678_i(y_i = data_help$Y, c_i = data_help$C,
k_i = data_help$K, x_i = data_help$X, l_i = data_help$L,
y_adj_bar = estimates[names(estimates) == "y_adj_bar"],
k_bar = mean(data_help$K),
alpha0 = estimates[names(estimates) == "alpha_0"],
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha2 = estimates[names(estimates) == "alpha_2"],
alpha3 = estimates[names(estimates) == "alpha_3"],
sigma1 = estimates[names(estimates) == "sigma_1"],
alpha4 = estimates[names(estimates) == "alpha_4"],
alphaXY = estimates[names(estimates) == "alpha_XY"],
sigma2sq = estimates[names(estimates) == "sigma_2_sq"]))
deriv_obj <- list(logL1_deriv = logL1_deriv, logL2_deriv = logL2_deriv)
}
return(deriv_obj)
}
#' @rdname score_and_hessian_matrix_functions
#' @export
scores <- function(derivobj = NULL) {
if (is.null(derivobj)) {
stop("derivobj has to be supplied.")
}
scores_out <- cbind(derivobj[[1]]$gradient[, 1:5],
derivobj[[2]]$gradient[, 6:8])
return(scores_out)
}
#' @rdname score_and_hessian_matrix_functions
#' @export
hessian <- function(derivobj = NULL) {
if (is.null(derivobj)) {
stop("derivobj has to be supplied.")
}
hessian_out <- derivobj[[1]]$hessian
hessian_out[, 6:8, ] <- derivobj[[2]]$hessian[, 6:8, ]
return(hessian_out)
}
Try the CIEE package in your browser
Any scripts or data that you put into this service are public.
CIEE documentation built on May 2, 2019, 6:39 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 deriv_obj scores hessian
Documented in deriv_obj hessian scores
#' Score and hessian matrix based on the estimating functions.
#'
#' Functions to compute the score and hessian matrices of the parameters
#' based on the estimating functions, under the GLM and AFT setting for
#' the analysis of a normally-distributed or censored time-to-event
#' primary outcome. The score and hessian matrices are further used in
#' the functions \code{\link{sandwich_se}}, \code{\link{ciee}} and
#' \code{\link{ciee_loop}} to obtain robust sandwich error estimates of the
#' parameter estimates of
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1^2, \alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1^2, \alpha4, \alphaXY, \sigma2^2}
#' under the GLM setting and
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1, \alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1, \alpha4, \alphaXY, \sigma2^2}
#' under the AFT setting.
#'
#' For the computation of the score and hessian matrices, first, the help function
#' \code{\link{deriv_obj}} is used. In a first step, the expression of all first
#' and second derivatives of the parameters is computed using the expressions of
#' \code{logL1} and \code{logL2} from the \code{\link{est_funct_expr}} as input.
#' Then, the numerical values of all first and second derivatives are obtained
#' for the observed data \code{Y}, \code{X}, \code{K}, \code{L} (and \code{C} under
#' the AFT setting) and point estimates (\code{estimates}) of the parameters,
#' for all observed individuals.
#'
#' Second, the functions \code{\link{scores}} and \code{\link{hessian}} are used
#' to extract the relevant score and hessian matrices with respect to \code{logL1}
#' and \code{logL2} from the output of \code{\link{deriv_obj}} and piece them together.
#' For further details, see the vignette.
#'
#' @param setting String with value \code{"GLM"} or \code{"AFT"} indicating
#' whether the matrices are computed under the GLM or AFT setting.
#' @param logL1 Expression of the function \code{logL1} generated by
#' the \code{\link{est_funct_expr}} function.
#' @param logL2 Expression of the function \code{logL2} generated by
#' the \code{\link{est_funct_expr}} function.
#' @param Y Numeric input vector for the primary outcome.
#' @param X Numeric input vector for the exposure variable.
#' @param K Numeric input vector for the intermediate outcome.
#' @param L Numeric input vector for the observed confounding factor.
#' @param C Numeric input vector for the censoring indicator under the AFT setting
#' (must be coded 0 = censored, 1 = uncensored).
#' @param estimates Numeric input vector with point estimates of the parameters
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1^2,}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1^2,}
#' \eqn{\alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha4, \alphaXY, \sigma2^2}
#' under the GLM setting and of
#' \eqn{\alpha_0, \alpha_1, \alpha_2, \alpha_3, \sigma_1, \alpha_4, \alpha_{XY}, \sigma_2^2}{\alpha0, \alpha1, \alpha2, \alpha3, \sigma1, \alpha4, \alphaXY, \sigma2^2}
#' under the AFT setting. Under the AFT setting,
#' \code{estimates} must also contain the mean of the estimated
#' true survival times \code{"y_adj_bar"}.
#' @param derivobj Output of the \code{\link{deriv_obj}} function used as input in
#' the \code{\link{scores}} and \code{\link{hessian}} functions.
#'
#' @return The \code{\link{deriv_obj}} function returns a list with
#' objects \code{logL1_deriv}, \code{logL2_deriv} which
#' contain the score and hessian matrices based on \code{logL1},
#' \code{logL2}, respectively.
#' @return The \code{\link{scores}} function returns the \eqn{(n \times 8)}{(n x 8)}
#' score matrix.
#' @return The \code{\link{hessian}} function returns the \eqn{(n \times 8 \times 8)}{(n x 8 x 8)}
#' hessian matrix.
#'
#' @examples
#'
#' # Generate data including Y, K, L, X under the GLM setting
#' dat <- generate_data(setting = "GLM")
#'
#' # Obtain estimating functions' expressions
#' estfunct <- est_funct_expr(setting = "GLM")
#'
#' # Obtain point estimates of the parameters
#' estimates <- get_estimates(setting = "GLM", Y = dat$Y, X = dat$X,
#' K = dat$K, L = dat$L)
#'
#' # Obtain matrices with all first and second derivatives
#' derivobj <- deriv_obj(setting = "GLM", logL1 = estfunct$logL1,
#' logL2 = estfunct$logL2, Y = dat$Y, X = dat$X,
#' K = dat$K, L = dat$L, estimates = estimates)
#' names(derivobj)
#' head(derivobj$logL1_deriv$gradient)
#'
#' # Obtain score and hessian matrices
#' scores(derivobj)
#' hessian(derivobj)
#'
#' @name score_and_hessian_matrix_functions
NULL
#' @rdname score_and_hessian_matrix_functions
#' @export
deriv_obj <- function(setting = "GLM", logL1 = NULL, logL2 = NULL, Y = NULL,
X = NULL, K = NULL, L = NULL, C = NULL,
estimates = NULL) {
if (is.null(setting) | is.null(logL1) | is.null(logL2) |
is.null(estimates)) {
stop("One or more arguments of the function are missing.")
}
if (is.null(Y) | is.null(X) | is.null(K) | is.null(L)) {
stop("Data of one or more variables are not supplied.")
}
if (setting == "AFT" & is.null(C)) {
stop("C has to be supplied for the AFT setting.")
}
n <- length(Y)
if (setting == "GLM") {
data_help <- data.frame(Y = Y, X = X, K = K, L = L)
data_help <- data_help[stats::complete.cases(data_help), ]
U12345_i <- stats::deriv(expr = logL1, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1sq", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "k_i", "x_i", "l_i", "alpha0",
"alpha1", "alpha2", "alpha3", "sigma1sq"), func = T,
hessian = T)
U678_i <- stats::deriv(expr = logL2, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1sq", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "k_i", "x_i", "y_bar", "k_bar",
"alpha1", "alpha4", "alphaXY", "sigma2sq"), func = T,
hessian = T)
logL1_deriv <- attributes(U12345_i(y_i = data_help$Y, k_i = data_help$K,
x_i = data_help$X, l_i = data_help$L,
alpha0 = estimates[names(estimates) == "alpha_0"],
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha2 = estimates[names(estimates) == "alpha_2"],
alpha3 = estimates[names(estimates) == "alpha_3"],
sigma1sq = estimates[names(estimates) == "sigma_1_sq"]))
logL2_deriv <- attributes(U678_i(y_i = data_help$Y, k_i = data_help$K,
x_i = data_help$X, y_bar = mean(data_help$Y),
k_bar = mean(data_help$K),
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha4 = estimates[names(estimates) == "alpha_4"],
alphaXY = estimates[names(estimates) == "alpha_XY"],
sigma2sq = estimates[names(estimates) == "sigma_2_sq"]))
deriv_obj <- list(logL1_deriv = logL1_deriv, logL2_deriv = logL2_deriv)
}
if (setting == "AFT") {
data_help <- data.frame(Y = Y, X = X, K = K, L = L, C = C)
data_help <- data_help[stats::complete.cases(data_help), ]
U12345_i <- stats::deriv(expr = logL1, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "c_i", "k_i", "x_i", "l_i",
"alpha0", "alpha1", "alpha2", "alpha3", "sigma1"),
func = T, hessian = T)
U678_i <- stats::deriv(expr = logL2, namevec = c("alpha0", "alpha1", "alpha2",
"alpha3", "sigma1", "alpha4", "alphaXY", "sigma2sq"),
function.arg = c("y_i", "c_i", "k_i", "x_i", "l_i",
"y_adj_bar", "k_bar", "alpha0", "alpha1", "alpha2",
"alpha3", "sigma1", "alpha4", "alphaXY", "sigma2sq"),
func = T, hessian = T)
logL1_deriv <- attributes(U12345_i(y_i = data_help$Y, c_i = data_help$C,
k_i = data_help$K, x_i = data_help$X, l_i = data_help$L,
alpha0 = estimates[names(estimates) == "alpha_0"],
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha2 = estimates[names(estimates) == "alpha_2"],
alpha3 = estimates[names(estimates) == "alpha_3"],
sigma1 = estimates[names(estimates) == "sigma_1"]))
logL2_deriv <- attributes(U678_i(y_i = data_help$Y, c_i = data_help$C,
k_i = data_help$K, x_i = data_help$X, l_i = data_help$L,
y_adj_bar = estimates[names(estimates) == "y_adj_bar"],
k_bar = mean(data_help$K),
alpha0 = estimates[names(estimates) == "alpha_0"],
alpha1 = estimates[names(estimates) == "alpha_1"],
alpha2 = estimates[names(estimates) == "alpha_2"],
alpha3 = estimates[names(estimates) == "alpha_3"],
sigma1 = estimates[names(estimates) == "sigma_1"],
alpha4 = estimates[names(estimates) == "alpha_4"],
alphaXY = estimates[names(estimates) == "alpha_XY"],
sigma2sq = estimates[names(estimates) == "sigma_2_sq"]))
deriv_obj <- list(logL1_deriv = logL1_deriv, logL2_deriv = logL2_deriv)
}
return(deriv_obj)
}
#' @rdname score_and_hessian_matrix_functions
#' @export
scores <- function(derivobj = NULL) {
if (is.null(derivobj)) {
stop("derivobj has to be supplied.")
}
scores_out <- cbind(derivobj[[1]]$gradient[, 1:5],
derivobj[[2]]$gradient[, 6:8])
return(scores_out)
}
#' @rdname score_and_hessian_matrix_functions
#' @export
hessian <- function(derivobj = NULL) {
if (is.null(derivobj)) {
stop("derivobj has to be supplied.")
}
hessian_out <- derivobj[[1]]$hessian
hessian_out[, 6:8, ] <- derivobj[[2]]$hessian[, 6:8, ]
return(hessian_out)
}
Try the CIEE package in your browser
Any scripts or data that you put into this service are public.
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.