R/samonCrossSummary.R

In samon: Sensitivity Analysis for Missing Data

# --------------------------------------------------------------------------------
# Summary function for samon differences in treatment for pairs of alpha
# --------------------------------------------------------------------------------
samonCrossSummary <- function(trt1,trt2, CIlevel=0.95) {

  Nalpha    <- trt1$Nalpha
  alphaList <- trt1$alphaList
  NSamples  <- trt1$NSamples
  n10       <- trt1$n0
  n20       <- trt2$n0
    
  TM1 <- trt1$TM[, c("alpha","IFEst","IFVar","meanjk","sejk","meanBoot","seBoot")]
  TM2 <- trt2$TM[, c("alpha","IFEst","IFVar","meanjk","sejk","meanBoot","seBoot")]
  TS1 <- trt1$TS[, c("alpha","Sample","IFEst","IFVar","meanjk","sejk","MIFEst")]
  TS2 <- trt2$TS[, c("alpha","Sample","IFEst","IFVar","meanjk","sejk","MIFEst")]

  colnames(TM1) <- c("alpha1","IFEst1","IFVar1","meanjk1","sejk1","meanBoot1","seBoot1")
  colnames(TM2) <- c("alpha2","IFEst2","IFVar2","meanjk2","sejk2","meanBoot2","seBoot2")
  colnames(TS1) <- c("alpha1","Sample","IFEst1","IFVar1","meanjk1","sejk1","MIFEst1")
  colnames(TS2) <- c("alpha2","Sample","IFEst2","IFVar2","meanjk2","sejk2","MIFEst2")

  for ( i1 in 1:Nalpha ) {
  for ( i2 in 1:Nalpha ) {
    alpha1 <- alphaList[i1]
    alpha2 <- alphaList[i2]

    TMA1 <- TM1[TM1[,"alpha1"] == alpha1,, drop=FALSE ]
    TMA2 <- TM2[TM2[,"alpha2"] == alpha2,, drop=FALSE ]
    TSA1 <- TS1[TS1[,"alpha1"] == alpha1,, drop=FALSE ]
    TSA2 <- TS2[TS2[,"alpha2"] == alpha2,, drop=FALSE ]
    
    TM                <- cbind(TMA1,TMA2)
    TM[,"Difference"] <- TM[,"IFEst2" ] - TM[,"IFEst1" ]
    TM[,"DIFSE"]      <- sqrt(TM[,"IFVar2"] + TM[,"IFVar1"])
    TM[,"DjkSE"]      <- sqrt(TM[,"sejk2"]^2 + TM[,"sejk1"]^2)
  
    TS <- merge(TSA1,TSA2,by=c("Sample"))
    TS[,"Difference"] <- TS[,"IFEst2" ] - TS[,"IFEst1" ]
    TS[,"DIFSE"]      <- sqrt(TS[,"IFVar2"] + TS[,"IFVar1"])
    TS[,"DjkSE"]      <- sqrt(TS[,"sejk2"]^2 + TS[,"sejk1"]^2)
    TS[,"TIF"]        <- (TS[,"Difference"] - TS[,"MIFEst2"] + TS[,"MIFEst1"]) / TS[,"DIFSE"]
    TS[,"TVals"]      <- (TS[,"Difference"] - TS[,"MIFEst2"] + TS[,"MIFEst1"]) / TS[,"DjkSE" ]
  
    TS    <- data.frame(TS)

    ## quantiles for q
  cicut  <- CIlevel
  cicut1 <- (1 - CIlevel)/2
  cicut2 <- 1 - cicut1
    myq <- function(x) {
      x1 <- as.vector(quantile(x, c(cicut1,cicut2)))
      x2 <- as.vector(quantile(abs(x), c(cicut)))
      return(c(x1,x2))                
    }
    ag1 <- matrix(myq( TS[,c("TVals")]), nrow=1 )
    ag2 <- matrix(myq( TS[,c("Difference")]), nrow=1)
    ag3 <- matrix(myq( TS[,c("TIF")  ]), nrow=1)

    colnames(ag1) <- c("tLow",   "tHigh",     "tSym")
    colnames(ag2) <- c("DIFLow", "DIFHigh",   "DIFSym")
    colnames(ag3) <- c("tIFLow", "tIFHigh",   "tIFSym")

    TM <- merge(TM,ag1)
    TM <- merge(TM,ag2)
    TM <- merge(TM,ag3)

    zcut <- qnorm( cicut2 )
    ## IF CI
    lb1 <- TM[,"Difference"] - zcut * TM[,"DIFSE"]
    ub1 <- TM[,"Difference"] + zcut * TM[,"DIFSE"]

    ## Bootstrap CI
    lb2 <- TM[,"Difference"] - zcut * TM[,"DIFSE" ]
    ub2 <- TM[,"Difference"] + zcut * TM[,"DIFSE" ]

    ## jk CI
    lb3 <- TM[,"Difference"] - zcut * TM[,"DjkSE" ]
    ub3 <- TM[,"Difference"] + zcut * TM[,"DjkSE" ]

    ## quantiles
    lb4 <- TM[,"DIFLow" ]
    ub4 <- TM[,"DIFHigh"]

    ## t using IF se
    lb5 <- TM[,"Difference"] - TM[,"tIFHigh"] * TM[,"DIFSE"]
    ub5 <- TM[,"Difference"] - TM[,"tIFLow" ] * TM[,"DIFSE"]

    ## t using jk se
    lb6 <- TM[,"Difference"] - TM[,"tHigh"] * TM[,"DjkSE"]
    ub6 <- TM[,"Difference"] - TM[,"tLow" ] * TM[,"DjkSE"]

    ## t using IF se symmetric
    lb7 <- TM[,"Difference"] - TM[,"tIFSym"] * TM[,"DIFSE"]
    ub7 <- TM[,"Difference"] + TM[,"tIFSym"] * TM[,"DIFSE"]
  
    ## t using jk se symmetric
    lb8 <- TM[,"Difference"] - TM[,"tIFSym"] * TM[,"DjkSE"]
    ub8 <- TM[,"Difference"] + TM[,"tIFSym"] * TM[,"DjkSE"]

    Difference <- TM[,"Difference"]
    CI <- cbind( alpha1, alpha2,Difference,lb1, ub1, lb2, ub2, lb3, ub3, lb4, ub4, lb5, ub5, lb6, ub6, lb7, ub7, lb8, ub8)

    if ( i1 == 1 & i2 == 1) {
       CIout <- CI
       TMout <- TM
    }
    else {
       CIout <- rbind(CIout,CI)
       TMout <- rbind(TMout,TM)
    }
}
}
 Ret <- list( TM = TMout, CI = CIout, Nalpha = Nalpha, alphaList = alphaList, CIlevel=CIlevel )
 return(Ret)
}