Merge pull request #24 from mbannick/develop

Develop

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: RobinCar
 Type: Package
 Title: Robust Estimation and Inference in Covariate-Adaptive Randomization
-Version: 0.2.0
+Version: 0.3.0
 Authors@R: c(
     person("Marlena", "Bannick", 
         email="[email protected]", 

NAMESPACE

@@ -3,22 +3,12 @@
 S3method(adjust,CovHR)
 S3method(adjust,CoxScore)
 S3method(adjust,GLMModel)
-S3method(adjust,LinModel)
 S3method(adjust,LogRank)
 S3method(adjust,SLModel)
 S3method(descript,GLMModelResult)
-S3method(descript,LinModelResult)
 S3method(jacobian,CUSTOM)
 S3method(jacobian,DIFF)
 S3method(jacobian,RATIO)
-S3method(linmod,ANCOVA)
-S3method(linmod,ANHECOVA)
-S3method(linmod,ANOVA)
-S3method(linmod,CUSTOM)
-S3method(lmlogic,ANCOVA)
-S3method(lmlogic,ANHECOVA)
-S3method(lmlogic,ANOVA)
-S3method(predictions,GLMModel)
 S3method(print,CalibrationResult)
 S3method(print,ContrastResult)
 S3method(print,GLMModelResult)
@@ -28,12 +18,6 @@ S3method(validate,RoboDataGLM)
 S3method(validate,RoboDataLinear)
 S3method(validate,RoboDataSL)
 S3method(validate,RoboDataTTE)
-S3method(vcov_car,GLMModel)
-S3method(vcov_car,LinModel)
-S3method(vcov_car,SLModel)
-S3method(vcov_sr,ANCOVA)
-S3method(vcov_sr,ANHECOVA)
-S3method(vcov_sr,ANOVA)
 export("%>%")
 export(car_pb)
 export(car_ps)
@@ -47,9 +31,7 @@ export(robincar_contrast)
 export(robincar_covhr)
 export(robincar_coxscore)
 export(robincar_glm)
-export(robincar_glm2)
 export(robincar_linear)
-export(robincar_linear2)
 export(robincar_logrank)
 export(robincar_tte)
 import(AIPW)
@@ -83,10 +65,10 @@ importFrom(dplyr,ungroup)
 importFrom(emulator,quad.tform)
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
-importFrom(stats,gaussian)
 importFrom(stats,model.matrix)
 importFrom(stats,rbinom)
 importFrom(stats,sd)
+importFrom(stats,setNames)
 importFrom(stats,uniroot)
 importFrom(stats,var)
 importFrom(survival,Surv)

NEWS.md

@@ -1,3 +1,20 @@
+# RobinCar 0.3.0
+
+## Features
+
+* Changed the name description of RobinCar.
+
+## Breaking Changes
+
+* Simplified user experience for `robincar_glm`: removed `robincar_glm`,
+and renamed `robincar_glm2` to `robincar_glm`. Same with `robincar_linear` and `robincar_linear2`. In effect,
+the older versions of `robincar_linear` and `robincar_glm` have been deprecated.
+
+## Bugfixes
+
+* Previously, the factor level order passed by the user for the treatment variable was ignored. We have fixed this,
+so all estimates will be presented in the order of treatment levels specified by user.
+
 # RobinCar 0.2.0
 
 ## Features

R/Data.R

@@ -116,6 +116,30 @@ validate.RoboDataTTE <- function(data, ref_arm, ...){
   return(list(newform, vars))
 }
 
+# Creates formula for robincar_linear
+# based on adj_method=c("ANOVA", "ANCOVA", "ANHECOVA").
+.create.formula <- function(adj_method, response_col, treat_col, covariate_cols){
+
+  if(adj_method == "ANOVA"){
+    form <- paste0(response_col, " ~ ", treat_col)
+  } else {
+
+    if(is.null(covariate_cols)) stop("Must provide covariates if
+                                     adjustment method is ANCOVA or ANHECOVA.")
+    covariates <- paste(covariate_cols, collapse=" + ")
+
+    if(adj_method == "ANCOVA"){
+      form <- paste0(response_col, " ~ ", treat_col, " + ", covariates)
+    } else if(adj_method == "ANHECOVA"){
+      form <- paste0(response_col, " ~ ", treat_col, " * (", covariates, ")")
+    } else {
+      stop("Unrecognized adjustment method.")
+    }
+
+  }
+  return(form)
+}
+
 .df.toclass <- function(df, classname, ...){
 
   data <- list()
@@ -146,7 +170,6 @@ validate.RoboDataTTE <- function(data, ref_arm, ...){
       # Include the covariates that will be needed for the formula
       # in a formula vector
       data[["formula"]] <- forms[[1]]
-      data[["formula_vars"]] <- df[forms[[2]]]
 
     } else if(grepl("col", att_name)){
 
@@ -176,13 +199,14 @@ validate.RoboDataTTE <- function(data, ref_arm, ...){
   data <- .df.toclass(df=df, classname=classname, ...)
 
   if(!is.null(data$treat)){
-    data$treat <- as.factor(as.vector(data$treat[[1]]))
+    data$treat <- data$treat[[1]]
+    if(!is.factor(data$treat)) data$treat <- as.factor(data$treat)
   }
   if(!is.null(data$response)){
-    data$response <- as.vector(data$response[[1]])
+    data$response <- data$response[[1]]
   }
   if(!is.null(data$event)){
-    data$event <- as.vector(data$event[[1]])
+    data$event <- data$event[[1]]
   }
   if(ncol(data$car_strata) == 0){
     data$car_strata <- NULL
@@ -197,10 +221,15 @@ validate.RoboDataTTE <- function(data, ref_arm, ...){
       data$car_strata[col] <- as.factor(data$car_strata[[col]])
     }
   }
-  if(ncol(data$covariate) == 0){
-    data$covariate <- NULL
+  if(!is.null(data$covariate)){
+    if(ncol(data$covariate) == 0){
+      data$covariate <- NULL
+    }
   }
 
+  # Save original data frame
+  data$df <- df
+
   # Add additional data attributes
   data$n <- nrow(df)
   data$k <- length(levels(data$treat))

R/adjust-glm.R

@@ -1,67 +1,81 @@
-predictions <- function(model, data, mod){
-  UseMethod("predictions", model)
-}
+# Fits a working model based on GLM and computes AIPW estimator
 
-#' @exportS3Method
-predictions.GLMModel <- function(model, data, mod, ...){
-  df <- data.frame(
-    treat=data$treat,
-    response=data$response
-  )
-  dmat <- get.dmat(data, model$adj_vars)
-  if(!is.null(dmat)){
-    df <- cbind(df, dmat)
+# This function allows us to implement the negative binomial model through MASS
+#' @importFrom MASS glm.nb
+fitmod <- function(family, ...){
+  useMASS <- FALSE
+  if(is.character(family)){
+    if(family == "nb"){
+      useMASS <- TRUE
+    }
   }
-  preds <- stats::predict(mod, newdata=df, type="response")
-  return(preds)
-}
-
-get.muhat <- function(model, data, mod){
-  set.treat <- function(a){
-    dat <- data
-    dat$treat <- rep(a, data$n)
-    dat$treat <- factor(dat$treat, levels=data$treat_levels)
-    return(dat)
+  if(useMASS){
+    mod <- MASS::glm.nb(...)
+    # copy over the 'model' attribute to the 'data'
+    # because it has a different name in glm.nb
+    mod$data <- mod$model
+  } else {
+    mod <- stats::glm(..., family=family)
   }
-  pred.treat <- function(dat) predictions.GLMModel(model, dat, mod)
-
-  datas <- lapply(data$treat_levels, set.treat)
-  preds <- lapply(datas, pred.treat)
-
-  pred_cols <- do.call(cbind, preds)
-  return(pred_cols)
+  return(mod)
 }
 
 # Perform GLM adjustment, based on the classes
 # of the model. Will perform adjustment based on the linear
 # model type of `model` and also do G-computation or AIPW
 # based on the second model type of `model`.
+#' @importFrom stats setNames
 #' @exportS3Method
 adjust.GLMModel <- function(model, data, ...){
 
-  # Get the generalized linear model and estimates by AIPW
-  glmod <- linmod(model, data, family=model$g_family, center=FALSE)
+  # 1. Set up data frame
+  df <- data$df
+  columns <- colnames(df)
+  columns[which(columns == data$treat_col)] <- "treat"
+  columns[which(columns == data$response_col)] <- "response"
+  df <- setNames(df, columns)
+
+  # 2. Fit GLM working model
+  glmod <- fitmod(
+    family=model$g_family,
+    stats::as.formula(data$formula),
+    data=df)
+
+  # 3. Predict potential outcome for each treatment group
+  # to compute \hat{\mu}_a for a = 1, ..., k
+  predict_a <- function(a){
 
-  # Get g-computation predictions
-  muhat <- get.muhat(model, data, glmod)
+    # Set up data frame
+    df_a <- df
+    df_a$treat <- rep(a, data$n)
+    df_a$treat <- factor(df_a$treat, levels=data$treat_levels)
+
+    # Create predictions
+    predictions <- stats::predict(glmod, newdata=df_a, type="response")
+
+    return(predictions)
+  }
+
+  preds <- lapply(data$treat_levels, predict_a)
+  muhat <- do.call(cbind, preds)
+
+  # 4. Save the g-computation estimate of \theta for diagnostic purposes
   g.estimate <- colMeans(muhat)
 
-  # Get mutilde from the GLM model, then estimate the treatment means by
-  # taking the mean over all of the potential outcomes
+  # 5. Use the \hat{\mu} to compute the adjusted \hat{\mu} for AIPW estimator
   glmod.adjusted <- get.mutilde(model, data, muhat)
   mutilde <- glmod.adjusted$mutilde
 
-  # Degree of freedom adjustment
-  df_adjust <- glmod.adjusted$df_adjust
-
+  # 6. Compute AIPW estimator of \theta
   estimate <- colMeans(mutilde)
 
-  # Compute the asymptotic variance
+  # 7. Compute the asymptotic variance
   asympt.variance <- vcov_car(model, data, glmod, mutilde)
-
+  df_adjust <- glmod.adjusted$df_adjust
   vcov_wt <- get.vcovHC(model$vcovHC, n=data$n, p=df_adjust)
   variance <- asympt.variance * vcov_wt / data$n
 
+  # 8. Format results
   result <- format_results(data$treat_levels, estimate, variance)
 
   return(