R/PAL.R

In ITRSelect: Variable Selection for Optimal Individualized Dynamic Treatment Regime

PAL <-
function(formula, data, subset, na.action, IC = c("BIC", "CIC", "VIC"),
             lambda.list = exp(seq(-3.5, 2, 0.1)), refit = TRUE, 
             control = PAL.control(...), model = TRUE, y = TRUE, a1 = TRUE, 
             x1 = TRUE, a2 = TRUE, x2 = TRUE, ...)
    {
        call <- match.call()
        if (missing(data))
            data <- environment(formula)
        # concordance information criterion is used by default	
        if (missing(IC))
            IC <- "CIC"
        # extract the model information	
        mf <- match.call(expand.dots = FALSE)
        m <- match(c("formula", "data", "subset", "na.action"), names(mf), 0L)
        mf <- mf[c(1L, m)]
        mf$drop.unused.levels <- TRUE
        oformula <- as.formula(formula)
        options(warn = -1)
        formula <- as.Formula(formula)	
        if (length(formula)[2L] < 4L){
            formula <- Formula(formula(formula, rhs = 1:2))	
        }
        else {
            formula <- Formula(formula(formula, rhs = 1:4))	
        }
        mf$formula <- formula
        mf[[1L]] <- as.name("model.frame")
        mf <- eval(mf, parent.frame())
        mt <- terms(formula, data = data)
        mtX1 <- delete.response(terms(formula, data = data, rhs = 1L))
        mtA1 <- delete.response(terms(formula, data = data, rhs = 2L))
        mtX2 <- delete.response(terms(formula, data = data, rhs = 3L))
        mtA2 <- delete.response(terms(formula, data = data, rhs = 4L))
        # obtain response
        Y <- model.response(mf, "numeric")
        # the first stage covariates
        X1 <- model.matrix(mtX1, mf)
		X1 <- as.matrix(X1[,-1])
        # the first stage treatment
        A1 <- model.matrix(mtA1, mf)
        A1 <- as.vector(A1[,-1])
        # the second stage covariates
        X2 <- model.matrix(mtX2, mf)
        # the second stage treatment
        A2 <- model.matrix(mtA2, mf)
        if (ncol(A2)>1){
            A2 <- as.vector(A2[,-1])
			X2 <- as.matrix(X2[,-1])
            # treatment is binary or not
            if (any((A2!=0) & (A2!=1)))
                stop("Treatment must be binary variable!")
        }
        else{
            A2 <- NULL
            X2 <- NULL
        }
        options(warn = 0)
        # treatment is binary or not
        if (any((A1!=0) & (A1!=1)))
            stop("Treatment must be binary variable!")
        if (length(Y) < 1)
            stop("empty model")
        result <- PAL.fit(y=Y, x1=X1, a1=A1, x2=X2, a2=A2, lambda.list=lambda.list, 
                          IC=IC, refit=refit, control=control)
        if (model)
            result$model <- mf
        if (y)
            result$y <- Y
        if (a1)
            result$a1 <- A1
        if (x1)
            result$x1 <- X1
        if (a2)
            result$a2 <- A2
        if (x2)
            result$x2 <- X2
		class(result) <- "PAL"
        return (result)	
    }