R/gem_fit.R

Defines functions gem_fit

Documented in gem_fit

#' GEM Fit
#' @name gem_fit
#' @import plyr
#' @import ggplot2
#' @useDynLib pirate
#' @importFrom Rcpp sourceCpp 
#' @importFrom stats lm rnorm
#' @importFrom RcppArmadillo fastLm
#' 
#' @description The main algorithm in \bold{pirate} package for calculating the coefficients of the linear combination of the covariates
#' to generate a GEM. This function can be applied to data sets with more than two treatment groups.
#' See more detail in \cite{E Petkova, T Tarpey, Z Su, and RT Ogden. Generated effect modifiers (GEMs) in randomized clinical trials. Biostatistics, (First published online: July 27, 2016). doi: 10.1093/biostatistics/kxw035.}
#'
#' @param dat Data frame with first column as the treatment index, second column as the outcome, and the 
#' remaining columns as the covariates design matrix. The elements of the treatment index take \eqn{K} distinct values, where \eqn{K} is the number of treatment groups. 
#' The outcome has to be a continuous variable.
#' @param method Choice of the criterion that the generated effect modifier optimizes. This is a string in 
#' \code{c("nu","de","F")}, which corresponde to the numerator, denominator and F-statistics
#' criteria respectively. The default method is the F-statistics method.

#' @return 
#' \enumerate{
#' 		\item \code{method} The criterion used to generate the GEM
#' 		\item \code{gemObject} Fitted result for the GEM model, see more in \bold{Details}
#' 		\item \code{p_value} The p-value for the interaction term in model \eqn{Y  = a + trt + Z + trt*Z + \epsilon}, where \eqn{Z} is the GEM
#' 		\item \code{Augmented_Data} The input data argumented with the GEM as the last column 
#' 		\item \code{effect.size} The effect size of the GEM if there are only two treatment groups
#' 		\item \code{plot} A scatter plot of Y versus the GEM with fitted lines and grouped by treatment
#' 		
#' }
#' 
#' @details \code{gemObject} is a list of three elements. The first element is the 
#' calculated weight \eqn{\alpha} for combining the predictors X. The second element contains the \eqn{K} vectors of coefficients \eqn{(\gamma_{j0},\gamma_{j1})} from  
#' \deqn{y_j = \gamma_{j0} + (X\alpha)\gamma_{j1} + \epsilon, j = 1,...,K, } for the \eqn{K} treatment groups respectively. The third element contains the \eqn{K}
#' vectors of coefficients from the unconstraint linear regression models \deqn{y_j = \beta_{j0} + X\beta_{j1} + \epsilon, j = 1,...,K, } for the \eqn{K} treatment groups respectively.
#'
#' @examples
#' #constructing the covariance matrix
#' co <- matrix(0.2, 10, 10)
#' diag(co) <- 1
#' dataEx <- data_generator1(d = 0.3, R2 = 0.5, v2 = 1, n = 300,
#'                         co = co, beta1 = rep(1,10),inter = c(0,0))
#' #fit the GEM
#' dat <- dataEx[[1]]
#' model_nu <- gem_fit(dat = dat, method = "nu")
#' model_de <- gem_fit(dat = dat, method = "de")
#' model_F <- gem_fit(dat = dat, method = "F")
#' @export
gem_fit <- function(dat, method = "F") 
{
       colnames(dat)[1:2] <- c("trt","y")
       k <- length(unique(dat$trt))
       dat_list <- dlply(dat,.(dat$trt),function(x){as.matrix(x[,-1])})

       gemObject <- gemCpp(dat_list, method)
       dat$Z <- as.matrix(dat[,c(-1,-2)]) %*% as.matrix(gemObject[[1]])
       
       mod <- fastLm(y~factor(trt)*Z,dat = dat)
       p_value <- summary(mod)[[2]][4,4]
       dat$trt <- factor(dat$trt)
       p <- ggplot(dat, aes_string(x = "Z", y = "y", group = "trt", color = "trt"))+
		       geom_point() + geom_smooth(method="lm",fullrange = T) 
        if ( k==2)
        {
            es <- effectSize(dat$y, dat$trt, dat$Z)
        }else
        {
            es <- NA
        }
	   results <- list("method"= method,
                       "gem"=gemObject,
                       "p_value" = p_value,
                       "Augmented_Data" = dat,
                       "effectSize" = es,
                       "plot" = p)
	return(results)	
	
}

Try the pirate package in your browser

Any scripts or data that you put into this service are public.

pirate documentation built on May 30, 2017, 2:22 a.m.