Nothing
#' Summarize a PSweight object
#'
#' \code{summary.PSweight} is used to summarize the results from \code{\link{PSweight}}.
#' The output contains the average causal effects defined by specific contrasts, as well as their
#' standard error estimates.
#'
#' @param object a PSweight object obtained from the \code{\link{PSweight}} function.
#' @param contrast a vector or matrix specifying the causal contrast of interest. The average causal effects will be
#' defined by such contrats. For multiple treatments, the contrast parameters are explained in Li and Li (2019)
#' for estimating general causal effects. Default is all pairwise contrasts between any two treatment groups.
#' @param type a character specifying the target estimand. The most commonly seen additive estimand is specified
#' by \code{type = "DIF"}, abbreviated for weighted difference-in-means. This is the usual pairwise average treatment
#' effects as those defined in Li, Morgan, and Zaslavsky (2018) and Li and Li (2019). For binary (or count outcomes), we also
#' allow two ratio estimands: causal relative risk (\code{type = "RR"}) and causal odds ratio (\code{type = "OR"}).
#' Estimates for these two ratio estimands will be reported on the log scale (log relative risk and log
#' odds ratio) to improve the approximate for asymptotic normality. With binary outcomes, \code{"DIF"} is the same
#' as the average causal risk difference. Default is "DIF" if left empty.
#' @param CI a logical argument indicates whether confidence interval should be calculated. Default is \code{CI = TRUE}.
#' @param ... further arguments passed to or from other methods.
#'
#' @details For the \code{contrast} argument, one specifies the contrast of interest and thus defines the target estimand
#' for comparing treatments. For example, if there are three treatment levels: A, B, and C, the contrast A-C
#' (i.e., E[Y(A)] - E[Y(C)]) can be specified by \code{c(1,0,-1)}. The contrasts of A-C and B-C can be
#' jointly specified by \code{rbind(c(1,0,-1), c(0,1,-1))}.
#'
#' For estimating the causal relative risk (\code{type = "RR"}), the contrast is specified at the log scale. For example,
#' the contrast A-C (specified by \code{c(1,0,-1)}) implies the estimation of log\{E[Y(A)]\} - log\{E[Y(C)]\}. For estimating the causal odds
#' ratio, the contrast is specified at the log odds scale. For example, the contrast A-C (specified by \code{c(1,0,-1)})
#' implies the estimation of log\{E[Y(A)]/E[1-Y(A)]\} - log\{E[Y(C)]/E[1-Y(C)]\}.
#'
#' The variance of the contrasts will be estimated by the delta method (if sandwich variance is used, or
#' \code{bootstrap = FALSE}), or nonparametric bootstrap (if \code{bootstrap = TRUE}). Details will be given in
#' Zhou et al. (2020+).
#'
#' The argument \code{type} takes one of three options: \code{"DIF"}, \code{"RR"}, or \code{"RR"}, with \code{"DIF"} as
#' the default option. Typically, \code{"RR"} is relavent for binary or count outcomes, and \code{"OR"} is relavent
#' only for binary outcomes. \code{"DIF"} applies to all types of outcomes.
#'
#' @return A list of following values:
#'
#' \describe{
#' \item{\code{ estimates}}{ a matrix of point estimates, standard errors, test statistics, 95% confidence intervals and p-values
#' for contrasts of interest.}
#'
#' \item{\code{ bootestimates}}{ a list of data frames containing estimated contrasts in each bootstrap replicate,
#' if bootstrap is used to estimate standard errors.}
#'
#' \item{\code{ contrast}}{ a table listing the specified contrasts of interest.}
#'
#' \item{\code{ group}}{ a table of treatment group labels corresponding to the output point estimates, provided in results
#' obtained from \code{\link{PSweight}}.}
#'
#' \item{\code{ trtgrp}}{ a character indicating the treatment group, or target population under ATT weights.}
#'
#' \item{\code{ type}}{ a character specifying the target estimand.}
#'
#' \item{\code{ CI}}{a logical indaicator of whether confidence interval should be reported.}
#' }
#'
#'
#' @references
#'
#' Li, F., Morgan, K. L., Zaslavsky, A. M. (2018).
#' Balancing covariates via propensity score weighting.
#' Journal of the American Statistical Association, 113(521), 390-400.
#'
#' Li, F., Li, F. (2019). Propensity score weighting for causal inference with multiple treatments.
#' The Annals of Applied Statistics, 13(4), 2389-2415.
#'
#'
#' @export
#'
#' @examples
#'
#' ## For examples, run: example(PSweight).
#'
#' @importFrom stats binomial coef cov formula glm lm model.matrix plogis poisson predict qnorm quantile sd
#' @importFrom utils capture.output combn
#' @importFrom graphics hist legend
#' @importFrom stats pnorm
summary.PSweight<-function(object,contrast=NULL,type='DIF',CI=TRUE,...){
muhat<-object$muhat
covmu<-object$covmu
muboot<-object$muboot
trtgrp<-object$trtgrp
group<-object$group
#groups
ngrp<-length(group)
#transform contrast into matrix
if(is.vector(contrast)){
contrast<-rbind(contrast)
}
#error message for wrong contrast
if(!is.null(contrast)){
if(dim(contrast)[2]!=ngrp){
cat('Contract length not equal to treatment groups, please check','\n')
cat('\n')
contrast<-NULL
}
cat('\n')
}
#if all contrast
if(is.null(contrast)){
ncst<-ngrp*(ngrp-1)/2
contrasttmp<-matrix(0,ngrp-1,ngrp)
contrasttmp[,1]<--1
contrasttmp[1:(ngrp-1),2:ngrp]<-diag(1,ngrp-1,ngrp-1)
contrast<-contrasttmp
if(ngrp>2){
for(i in 1:(ngrp-2)){
ntmp<-nrow(contrasttmp)-1
contrasttmp<-cbind(0,contrasttmp)[1:ntmp,1:ngrp]
contrast<-rbind(contrast,contrasttmp)
}
}
}
rownames(contrast)<-paste('Contrast',1:nrow(contrast))
colnames(contrast)<-group
#define the bootstrap p value function
pvalboot <- function(x,est) {
esttmp<-abs(est)
x<-abs(x-mean(x))
return(mean(x>esttmp))
}
#use bootstrap or not
closedform<-is.null(object$muboot)
if(closedform){
if(type=='DIF'){
est.h<-c(contrast%*%muhat)
se.h<-sqrt(diag(contrast%*%covmu%*%t(contrast)))
zval.h<-est.h/se.h
lcl<-est.h-qnorm(0.975)*se.h
ucl<-est.h+qnorm(0.975)*se.h
p.value<-2*(1-pnorm(abs(est.h/se.h)))
}else if(type=='RR'){
tranmuhat<-log(muhat)
tranest<-c(contrast%*%tranmuhat)
trancovmu<-diag(1/muhat)%*%covmu%*%diag(1/muhat)
transe<-sqrt(diag(contrast%*%trancovmu%*%t(contrast)))
tranlcl<-tranest-qnorm(0.975)*transe
tranucl<-tranest+qnorm(0.975)*transe
est.h<-(tranest)
se.h<-transe
zval.h<-est.h/se.h
p.value<-2*(1-pnorm(abs(est.h/se.h)))
lcl<-(tranlcl)
ucl<-(tranucl)
}else if(type=='OR'){
tranmuhat<-log(muhat/(1-muhat))
tranest<-c(contrast%*%tranmuhat)
trancovmu<-diag(1/(muhat*(1-muhat)))%*%covmu%*%diag(1/(muhat*(1-muhat)))
transe<-sqrt(diag(contrast%*%trancovmu%*%t(contrast)))
tranlcl<-tranest-qnorm(0.975)*transe
tranucl<-tranest+qnorm(0.975)*transe
est.h<-(tranest)
se.h<-transe
zval.h<-est.h/se.h
p.value<-2*(1-pnorm(abs(est.h/se.h)))
lcl<-(tranlcl)
ucl<-(tranucl)
}else{
cat('type not found')
cat('\n')
}
}else{
if(type=='DIF'){
samp<-muboot%*%t(contrast)
est.h<-c(contrast%*%muhat)
se.h<-sqrt(diag(cov(samp)))
zval.h<-est.h/se.h
lcl<-apply(samp,2,function(x) quantile(x,0.025))
ucl<-apply(samp,2,function(x) quantile(x,0.975))
p.value<-c()
dimcon<-dim(samp)[2]
for(j in 1:dimcon){
esttmp<-c(est.h)[j]
p.value<-c(p.value,pvalboot(samp[,j],esttmp))
}
}else if(type=='RR'){
samp<-log(muboot)%*%t(contrast)
tranmuhat<-log(muhat)
tranest<-c(contrast%*%tranmuhat)
est.h<-tranest
se.h<-sqrt(diag(cov(samp)))
zval.h<-est.h/se.h
lcl<-apply(samp,2,function(x) quantile(x,0.025))
ucl<-apply(samp,2,function(x) quantile(x,0.975))
p.value<-c()
dimcon<-dim(samp)[2]
for(j in 1:dimcon){
esttmp<-c(est.h)[j]
p.value<-c(p.value,pvalboot(samp[,j],esttmp))
}
}else if(type=='OR'){
samp<-log(muboot/(1-muboot))%*%t(contrast)
tranmuhat<-log(muhat/(1-muhat))
tranest<-c(contrast%*%tranmuhat)
est.h<-tranest
se.h<-sqrt(diag(cov(samp)))
zval.h<-est.h/se.h
lcl<-apply(samp,2,function(x) quantile(x,0.025))
ucl<-apply(samp,2,function(x) quantile(x,0.975))
p.value<-c()
dimcon<-dim(samp)[2]
for(j in 1:dimcon){
esttmp<-c(est.h)[j]
p.value<-c(p.value,pvalboot(samp[,j],esttmp))
}
}else{
cat('type not found')
cat('\n')
}
}
estimates<-cbind(est.h,se.h,zval.h,lcl,ucl,p.value)
if(is.vector(estimates)){
t(estimates)
}
colnames(estimates)<-c("Estimate","Std.Error","z value", "lwr","upr","Pr(>|z|)")
rownames(estimates)<-rownames(contrast)
if(is.null(object$muboot)){
bootestimates<-NULL
}else{
if(is.vector(samp)){
samp<-t(samp)
}
bootestimates<-samp
colnames(bootestimates)<-rownames(contrast)
rownames(bootestimates)<-NULL
}
out<-list(estimates=estimates,bootestimates=bootestimates,contrast=contrast,group=group,trtgrp=trtgrp,type=type,CI=CI)
class(out)<-'PSweightsum'
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.