R/sandwich.r

In spass: Study Planning and Adaptation of Sample Size

#' @title Calculate the robust covariance estimator for GEE given an
#'
#' @description \code{sandwich} calculates the covariance structure between timepoints given matrices \code{yCov}, \code{D},\code{V} and \code{correctionmatrix}. This is done to be able to account for missingness in the Data.
#'
#' @param yCov              \eqn{yCov} matrix containing either an estimation for the covariance between timepoints or an empirical covariance matrix itsel. see 'Details'.
#' @param D                 \eqn{D} denotes the mean matrix of all entries of \eqn{\Delta\mu_i/\delta\beta}, where this is the average over all \code{i} patiens. see 'Details'.
#' @param V                 \eqn{V} denotes the working covariance matrix.  see 'Details'.
#' @param missing           vector which denotes the probability to experience a dropout at each timepoint. If \code{missingtype} is \code{"none"} then all entries are 0.
#' @param missingtype       String which describes the type of missingness occuring in the data. \code{none} if no missingnes occured, \code{"monotone"} if missing was monotone and \code{"intermittened"} if the missingness was independent across all timepoints.
#' @param correctionmatrix  As of this version this matrix is needed to correct some calculations. see 'Details' to see for more details and how to correctly select matrices.
#'
#' @details \code{yCov} is either empirical or the estimated covariance-matrix between timepoints which is needed to calculate the sandwich estimator. This matrix can either be generated by estimating the empirical covariance matrix using existing data or by using function \code{gen_cov_cor} to calculate a estimation for the covariance.
#'
#' \code{D} denotes the estimation of \eqn{n^-1* \sum_i^N \Delta\mu_i/\delta\beta}, which means that \eqn{D=E(D_i)}. As of yet this has the unfortunate side effect that E(D_i%*%D_i) is not equal to E(D_i)%*%E(D_i). Entries that differ between those two formulas are equated using \code{correctionsmatrix}.
#'
#'
#' @return \code{sandwich} returns the robust covariance estimator of regression coefficients which are implicitly defined by \code{D}.
#'
#' @source \code{sandwich} computes the asymptotic sandwich covariance estimator and uses code contributed by Roland Gerard Gera.
#'
#' @references Liang Kung-Yee, Zeger Scott L. (1986); Jung Sin-Ho, Ahn Chul (2003); Wachtlin Daniel Kieser Meinhard (2013)
#'
#' @examples
#' #Let's assume we wish to calculate the robust variance estimator for  equation
#' #\eqn{y_{it}=\beta_0+\beta_1*I_{treat}+\beta_2*t+\beta_3*I _{treat}*t+\epsilon_{it}}.
#' #Furthermore we use the identitiy matrix as the working covariance matrix.
#' #The chance to get treatment is 60 percent and the observed timerange ranges from 0:5.
#'
#'   ycov = gen_cov_cor(var = 3,rho = 0.25,theta = 1,Time = 0:5,cov = TRUE)
#'   D = matrix(c(1,0.6,0,0,
#'                1,0.6,1,0.6,
#'                1,0.6,2,1.2,
#'                1,0.6,3,1.8,
#'                1,0.6,4,2.4,
#'                1,0.6,5,3.0),nrow=4)
#'
#'  D=t(D)
#'  V=diag(1,length(0:5))
#'  #We correct entries where E(D_i %*% D_i) is unequal to E(D_i)%*%E(D_i) (D %*% D).
#'  correctionmatrix=matrix(c(1,1,1,1,1,1/0.6,1,1/0.6,1,1,1,1,1,1/0.6,1,1/0.6),nrow=4)
#'  missingtype = "none"
#'
#'  robust=sandwich(yCov=ycov,D=D,V=V,missingtype=missingtype,correctionmatrix=correctionmatrix)
#'  robust
#'
#' @import MASS
#' @export

sandwich <- function(yCov , D, V, correctionmatrix,missing=rep(0,dim(yCov)[[2]]), missingtype=c("none","monotone","intermittened") ){
  #  require(MASS)

  if(missingtype=="none"){
    remain_i = rep(1,dim(yCov)[[2]])
    remain_ij=remain_i %*% t(remain_i) #p_{i,j}= p_i*p_j
  }
  if(missingtype=="intermittened"){
    # calculate percentage of remaining subjects per timepoint
    remain_i = 1-missing
    # calculate probabilety to remain in study at time i AND j
    remain_ij=remain_i %*% t(remain_i) #p_{i,j}= p_i*p_j
    diag(remain_ij)=remain_i # bzw p_{i,i}=\p_i
  }
  if(missingtype=="monotone"){
    # calculate percentage of remaining subjects per timepoint
    remain_i = 1-missing
    # calculate probabilety to remain in study at time i AND j
    remain_ij=diag(rep(0,dim(yCov)[[1]]))
    for (i in 1:length(missing)){
      for (j in 1:length(missing)){
        # alternativly calculate min(remain_i[i],remain_i[j])]
        remain_ij[i,j]=remain_i[max(i,j)]
      }
    }
  }

  Bread=t(D)%*%diag(remain_i)%*%ginv(V)%*%D*correctionmatrix
  invbread=ginv(Bread)

  Meat=(t(D)%*%ginv(V)%*%(remain_ij*yCov)%*%ginv(V)%*%D)*correctionmatrix

  robust=invbread%*%Meat%*%invbread
  return(robust)
}