Monte.Carlo.se (version 0.1.1)

mc.se.matrix: Standard Errors for Summaries Based on 2 or More Vectors of Monte Carlo Output

Description

mc.se.matrix --- jackknife and bootstrap SEs for statistics based on k correlated samples

Usage

mc.se.matrix(x, xcol, B = 0, seed = NULL, summary.f, ...)

Value

Returns data frame of summary.f, Monte Carlo (MC) SE of summary.f, MC sample size N, method (jackknife or bootstrap), B and seed if bootstrap is used

Arguments

x

N by k matrix or data frame of MC output, k>1

xcol

vector specifying columns of x to use

B

B=0 (default) means use jackknife, B>0 means use bootstrap with B resamples, If B>0, then a seed must be given to start the bootstrap resampling

seed

seed=NULL (default) used with jackknife, otherwise needs positive integer

summary.f

summary function with arguments (index,xdata), x goes into xdata

...

Additional arguments to be passed

Author

Dennis Boos, Kevin Matthew, Jason Osborne

Details

Suppose that k columns of Monte Carlo output are produced from N independent Monte Carlo samples, and summary statistics involving 2 or more columns are to be reported in tables. mc.se.matrix gives Monte Carlo standard errors (SEs) for these summary statistics (e.g., the ratio of 2 column means or variances). The vignette vignette("Example2", package = "Monte.Carlo.se") is a detailed account of using mc.se.vector.

References

Boos, D. D., and Osborne, J. A. (2015), "Assessing Variability of Complex Descriptive Statistics in Monte Carlo Studies using Resampling Methods," International Statistical Review, 25, 775-792.

Examples

# \donttest{
# First create MC output used in Table 9.1, p. 367, of Boos and Stefanski (2013).
# norm15 holds 10,000 sample means, 20% trimmed means, and medians
# for normal samples of size 15

N <- 10000
set.seed(346)                   # sets the random number seed
z <- matrix(rnorm(N*15),nrow=N) # N rows of N(0,1) samples, n=15
out.m.15 <- apply(z,1,mean)     # mean for each sample
out.t20.15 <- apply(z,1,mean,trim=0.20)   # 20% trimmed mean for each sample
out.med.15 <- apply(z,1,median)           # median for each sample

# Save all 1000 blocks of 3 estimators in a data frame
norm15 <- data.frame(mean=out.m.15,trim20=out.t20.15,median=out.med.15)

# Pearson correlation based on 2 vectors of MC output
# Note that the 2 vectors are in xdata, index is for the rows of xdata
corr<-function(index,xdata){cor(xdata[index,1],xdata[index,2])}

# Compute jackknife SE for summary.f=corr

mc.se.matrix(norm15,xcol=c(1,2),summary.f=corr)
#      summary          se     N    method
#  1 0.9367602 0.001256079 10000 Jackknife

# Compute bootstrap SE for summary.f=corr

mc.se.matrix(norm15,xcol=c(1,2),summary.f=corr,B=1000,seed=3928)
#      summary          se     N    method    B seed
#  1 0.9367602 0.001287065 10000 Bootstrap 1000 3928

# Rerun with B=5000

mc.se.matrix(norm15,xcol=c(1,2),summary.f=corr,B=5000,seed=3928)
#      summary          se     N    method    B seed
#  1 0.9367602 0.001266177 10000 Bootstrap 5000 3928

# Compute jackknife SE for summary.f=ratio.var
# = ratio of variances of the two columns
# A ratio of 2 variances facilitates comparison of the variances

ratio.var <- function(index,xdata)
{var(xdata[index,1])/var(xdata[index,2])}

# ratio of column 1 variance to column 2 variance

mc.se.matrix(norm15,xcol=c(1,2),summary.f=ratio.var)
#     summary          se     N    method
# 1 0.8895367 0.006263652 10000 Jackknife
# Coupled with SE=0.006, the ratio=0.89 shows the second variance is larger than the fitst

# ratio of column 2 variance to column 1 variance
# Same conclusion as for the previous ratio

mc.se.matrix(norm15,xcol=c(2,1),summary.f=ratio.var)
#    summary          se     N    method
# 1 1.124181 0.007915667 10000 Jackknife
# }

# Function Code

mc.se.matrix <- function(x,xcol,B=0,seed=NULL,summary.f,...){
 # x is an n by k matrix or data frame of MC output, k>1
 # xcol is a vector specifying columns of x to use
 # summary.f is the summary function with arguments (index,xdata)
 # ... is for additional arguments to summary.f
 # B=0 means use jackkife, B>0 means use bootstrap with B resamples
 # If B>0, then a seed must be given to start the bootstrap resampling
 N=nrow(x)
 x=x[,xcol]   # use columns selected
 if(B>0){
   if(is.null(seed))stop('If B>0, then seed for the bootstrap must be given in the call')
   se=boot.se(1:N,B=B,theta = summary.f,xdata=x,...)}
 else{se=jack.se(1:N,theta = summary.f,xdata=x,...)}
 summ = summary.f(1:N,xdata=x,...)
 if(B>0) {out=data.frame(summary=summ,se,N,method="Bootstrap",B,seed)}
 else{out=data.frame(summary=summ,se,N,method="Jackknife")}
 return(out)}