[R] Bootstrap CIs for weighted means of paired differences

i.petzev i.petzev at gmail.com
Thu Nov 20 11:23:30 CET 2014


Hi David,

sorry, I was not clear. The difference comes from defining or not defining “w” in the boot() function. The results with your function and your approach are thus:

set.seed(1111)
x <- rnorm(50)
y <- rnorm(50)
weights <- runif(50)
weights <- weights / sum(weights)
dataset <- cbind(x,y,weights)

vw_m_diff <- function(dataset,w) {
  differences <- dataset[w,1]-dataset[w,2]
  weights <- dataset[w, "weights"]
  return(weighted.mean(x=differences, w=weights))
}
res_boot <- boot(dataset, statistic=vw_m_diff, R = 1000, w=dataset[,3])
boot.ci(res_boot)

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res_boot)

Intervals : 
Level      Normal              Basic         
95%   (-0.5657,  0.4962 )   (-0.5713,  0.5062 )  

Level     Percentile            BCa          
95%   (-0.6527,  0.4249 )   (-0.5579,  0.5023 )  
Calculations and Intervals on Original Scale

********************************************************************************************************************

However, without defining “w” in the bootstrap function, i.e., running an ordinary and not a weighted bootstrap, the results are:

res_boot <- boot(dataset, statistic=vw_m_diff, R = 1000)
boot.ci(res_boot)

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res_boot)

Intervals : 
Level      Normal              Basic         
95%   (-0.6265,  0.4966 )   (-0.6125,  0.5249 )  

Level     Percentile            BCa          
95%   (-0.6714,  0.4661 )   (-0.6747,  0.4559 )  
Calculations and Intervals on Original Scale

On 19 Nov 2014, at 17:49, David Winsemius <dwinsemius at comcast.net> wrote:

>>> vw_m_diff <- function(dataset,w) {
>>>     differences <- dataset[w,1]-dataset[w,2]
>>>    weights <- dataset[w, "weights"]
>>>    return(weighted.mean(x=differences, w=weights))
>>>  }


	[[alternative HTML version deleted]]



More information about the R-help mailing list