Nothing
get_design_matrix <- function(
formula, data, response = TRUE, predictors = TRUE,
Terms = NULL, has_fixed_intercept = NULL, xlev = NULL, fixed = NULL,
random = NULL){
if(is.null(fixed) && is.null(random))
return(.get_design_matrix_one_frm (
formula = formula, data = data, response = response,
predictors = predictors, Terms = Terms,
has_fixed_intercept = has_fixed_intercept, xlev = xlev))
.get_design_matrix_two_frms(
data = data, response = response, predictors = predictors, Terms = Terms,
xlev = xlev, fixed = fixed, random = random)
}
# find the design matrix and returns the left hand side and right hand side of
# objects in the formula
.get_design_matrix_one_frm <- function(
formula, data, response = TRUE, predictors = TRUE,
Terms = NULL, has_fixed_intercept = NULL, xlev = NULL){
if(!response && !predictors)
stop(sQuote("get_design_matrix"),
" is called where neither respone or predictors is requested")
if(!all(
is.null(Terms) == c(is.null(has_fixed_intercept), is.null(xlev))))
stop(sQuote("Terms"), ", ", sQuote("has_fixed_intercept"), " and ",
sQuote("xlev"), " must either all be NULL or not")
if(is_first_call <- is.null(Terms)){
#####
# See if there is a fixed intercept. If so, update the formula
tt <- terms(formula, data = data, specials = "ddFixed_intercept")
fixed_inter_terms <- attr(tt, "specials")[["ddFixed_intercept"]]
if(length(fixed_inter_terms) > 1)
stop("more than one ", sQuote("ddFixed_intercept"), " term")
has_fixed_intercept <- length(fixed_inter_terms) > 0
formula_used <- if(has_fixed_intercept){
attr(tt, "intercept") <- 1
tt <- drop.terms(tt, fixed_inter_terms - attr(tt, "response"),
keep.response = TRUE)
attributes(tt) <- attributes(formula)
tt
} else
formula
}
do_keep_random_intercept <- .keep_random_intercept(formula, data = data)
#####
# Get model.frame, model.matrix and outcome. The two latter if needed
if(is_first_call){
Call <- match.call()
indx <- match(c("formula", "data"), names(Call), nomatch = 0)
temp <- Call[c(1, indx)]
temp[[1]] <- as.name("model.frame")
formula_used <- terms(formula_used, data = data, specials = "ddFixed")
temp[[2]] <- formula_used
if(!response)
temp$formula <- delete.response(temp$formula)
if(!predictors)
temp$formula <- update(temp$formula, . ~ -1)
mf <- eval(temp, parent.frame())
} else {
if(!response)
Terms <- delete.response(Terms)
if(!predictors)
Terms <- update(Terms, . ~ -1)
mf <- model.frame(Terms, data, xlev = xlev)
formula_used <- Terms
}
Y <- if(response) model.response(mf) else NULL
if(!predictors){
Terms <- NULL
X <- NULL
fixed_terms <- NULL
xlev <- NULL
} else {
Terms <- terms(mf)
X <- model.matrix(Terms, mf)
xlev <- .getXlevels(Terms, mf)
fixed_terms_indicies <- attr(Terms, "specials")[["ddFixed"]]
if(is.null(fixed_terms_indicies)){
fixed_terms_indicies <- c()
} else {
fixed_terms_indicies <- which(
attr(X, "assign") %in%
(fixed_terms_indicies - attr(Terms, "response")))
}
fixed_terms <- X[, fixed_terms_indicies, drop = FALSE]
if(length(fixed_terms_indicies) > 0)
X <- X[, -fixed_terms_indicies, drop = FALSE]
if(has_fixed_intercept){
fixed_terms <- cbind(`(Intercept)` = 1, fixed_terms)
if(!do_keep_random_intercept)
X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
}
}
list(X = X, fixed_terms = fixed_terms, Y = .set_Y_to_counting(Y),
formula_used = formula_used,
terms = Terms, has_fixed_intercept = has_fixed_intercept, xlev = xlev)
}
.keep_random_intercept <- function(formula, data){
tt <- terms(formula, specials = "ddFixed_intercept", data = data)
fixed_inter_terms <- attr(tt, "specials")[["ddFixed_intercept"]]
if(length(fixed_inter_terms) == 0)
return(FALSE)
if(length(fixed_inter_terms) > 1)
stop("more than one ", sQuote("ddFixed_intercept"), " term")
fixed_expr <- attr(tt, "variables")[[fixed_inter_terms + 1L]]
attr(eval(fixed_expr), "random_intercept")
}
.get_design_matrix_two_frms <- function(
data, response = TRUE, predictors = TRUE, Terms = NULL, xlev = NULL, fixed = NULL,
random = NULL){
if(!is.null(Terms) || !is.null(xlev))
stop("method not implemented")
stopifnot(inherits(fixed, "formula"), inherits(random, "formula"))
. <- function(formula, data, ...){
tt <- terms(formula, data, specials = c("ddFixed", "ddFixed_intercept"))
if(any(sapply(attr(tt, "specials"), length) > 0))
stop("Do not use ", sQuote("ddFixed"), " or ",
sQuote("ddFixed_intercept"),
" when you use ", sQuote("random"), " and ", sQuote("fixed"))
cl <- match.call()
cl[[1L]] <- quote(model.frame)
eval(cl, parent.frame())
}
if(!predictors)
fixed <- update(fixed, . ~ 1)
mf <- .(fixed, data)
Y <- if(response)
model.response(mf) else NULL
if(predictors){
Terms_fixed <- terms(mf)
xlev_fixed <- .getXlevels(Terms_fixed, mf)
fixed_terms <- model.matrix(Terms_fixed, mf)
mf <- model.frame(random, data)
Terms_random <- terms(mf)
xlev_random <- .getXlevels(Terms_random, mf)
X <- model.matrix(Terms_random, mf)
} else
X <- fixed_terms <- Terms_fixed <- Terms_random <- xlev_fixed <-
xlev_random <- NULL
list(X = X, fixed_terms = fixed_terms, Y = .set_Y_to_counting(Y),
terms = list(fixed = Terms_fixed, random = Terms_random),
xlev = list(fixed = xlev_fixed , random = xlev_random))
}
.set_Y_to_counting <- function(Y){
# Change outcome if formula it was Surv(stop, event) ~
if(!is.null(Y) && attr(Y, "type") == "right"){
Y <- cbind(0, Y)
dimnames(Y)[[2]] <- c("start", "stop", "status")
attr(Y, "type") <- "counting"
}
Y
}
#' Auxiliary Functions for Fixed Effects
#' @description
#' Functions used in formula of \code{\link{ddhazard}} for time-invariant effects. \code{ddFixed_intercept} is only used for the intercept.
#'
#' @param object expression that would be used in formula. E.g. \code{x} or \code{poly(x, degree = 3)}.
#'
#' @return
#' Returns the passed object.
#'
#' @examples
#' # we can get a time-invariant effect of `x1` by
#' set.seed(1)
#' dat <- data.frame(stop = 1:20, event = rep(c(TRUE, FALSE), 10L), x1 = rnorm(20))
#' ddhazard(Surv(stop, event) ~ ddFixed(x1), dat,
#' Q_0 = diag(1), by = 1, Q = diag(1))
#'
#' # all of the calls below will yield the same result with a time-invariant
#' # intercept:
#' ddhazard(Surv(stop, event) ~ ddFixed_intercept() + x1, dat,
#' Q_0 = diag(1), by = 1, Q = diag(1))
#' ddhazard(Surv(stop, event) ~ -1 + ddFixed_intercept() + x1, dat,
#' Q_0 = diag(1), by = 1, Q = diag(1))
#' @export
ddFixed <- function(object){
if(all(object == 1)){
# check that the function is not called as part of predict where we may
# just have ones
funcs <- lapply(sys.calls(), "[[", 1)
if(!any(sapply(funcs, function(x)
any(sapply(x, function(z) is.name(z) && z == as.name("predict"))))))
stop("All elements in call to ", sQuote("ddFixed"), " is one.",
" This is potentially a depreciated way of fixing the intercept.",
" See ", sQuote("?ddFixed"))
}
object
}
#' @rdname ddFixed
#'
#' @param random_intercept \code{TRUE} if a zero mean time-varying process
#' should be included at as an additional term. Only relevant in stationary
#' models. See the \code{type} argument in \code{\link{PF_EM}}.
#'
#' @export
ddFixed_intercept <- function(random_intercept = FALSE){
if(!is.logical(random_intercept) || length(random_intercept) != 1L)
stop(sQuote("random_intercept"), " needs to be a scalar logical")
structure(list(), random_intercept = random_intercept)
}
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.