Bootstrapping RUM and RRM models with no variable splits.

rm(list = ls(all = TRUE))
library(ggplot2)
### RUM general results
load('RUM_wtp_boot.RData')
rum_wtp <- sapply(1:length(RUM_wtp_boot), function(x) RUM_wtp_boot[[x]][['wtp']])

### RRM General Results
load('RRM_wtp_boot.RData')
wtp_opt_out = as.vector(sapply(RRM_wtp_boot, mean))

Data reshaping: Take the mean of each vector of WTP from the RRM model

mean_wtp_opt_out_rrm <- data.frame(value=do.call(c, lapply(wtp_opt_out, function(x) mean(x))))
mean_wtp_opt_out_rrm$variable <- 'RRM'
rum_wtp <- data.frame(value=rum_wtp)
rum_wtp$variable <- 'RUM'

# Prepare the data for the plot
mdata <- data.frame(rbind(mean_wtp_opt_out_rrm, rum_wtp))
names(mdata) <- c('Mean_WTP', 'Model')
mdata$Model <- as.factor(mdata$Model)
p <- ggplot(mdata, aes(x=Model, y=Mean_WTP, fill=Model))
p <- p + geom_boxplot()
p <- p + theme_bw()
p <- p + scale_fill_manual(values=c('grey40', 'grey70'))
p <- p + labs(list(title="Comparison between RRM and RUM model's average willingness to pay",
                   y = 'Average willingness to pay'))
print(suppressWarnings(p))

p <- ggplot(mdata, aes(x=Mean_WTP, fill=Model, color=Model))
p <- p + geom_histogram(aes(y = ..density..), alpha=.5) + geom_density(alpha=0.3)
p <- p + scale_fill_manual(values=c('grey40', 'grey70'))
p <- p + theme_bw()
p <- p + labs(list(title="Comparison between RRM and RUM model's average willingness to pay",
                   x = 'Bootstrapped Willingness to pay', y='Density'))
print(suppressWarnings(p))

distribution comparison

wtp_rum <- sort(mdata[mdata$Model == 'RUM', 'Mean_WTP'])
wtp_rrm <- sort(mdata[mdata$Model == 'RRM', 'Mean_WTP'])

` KS and Wilcox tests

#library(plotrix)
ks.test(wtp_rrm, wtp_rum)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm and wtp_rum
## D = 0.48, p-value = 1.387e-05
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm, wtp_rum)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm and wtp_rum
## W = 1827, p-value = 7.059e-05
## alternative hypothesis: true location shift is not equal to 0
std.error <- function(x) sd(x)/sqrt(length(x))

# summary functions used below
sumarizar <- function(x){
  ls = list(c(mean(x), std.error(x), sd(x)), quantile(x, probs=c(0, 0.05, 0.25, 0.5, 0.75, 0.95, 1)))
  names(ls)[[1]] <- 'Mean - SE - SD'
  names(ls)[[2]] <- 'Quantile'
  
  conf_interval <- c(mean(x) - 1.96 * std.error(x), mean(x) + 1.96 * std.error(x))
  ls$conf_inf <- conf_interval
  names(ls)[[3]] <- 'Confidence interval'
  return(ls)
}

mean_density <- function(wtp_rum_i, avgwtp){
  ns <- length(wtp_rum_i)
  density <- rep(NA, ns)
  for(i in 1:ns){
    density[i] <- mean( (wtp_rum_i[i] > avgwtp))
  }
  mean(density)
}

Rest of the summary results

sumarizar(wtp_rum)
## $`Mean - SE - SD`
## [1] 104.065352   4.307416  30.458034
## 
## $Quantile
##        0%        5%       25%       50%       75%       95%      100% 
##  22.63914  50.97024  86.28236 110.62524 123.24312 146.63545 158.54812 
## 
## $`Confidence interval`
## [1]  95.62282 112.50789
sumarizar(wtp_rrm)
## $`Mean - SE - SD`
## [1] 131.558199   5.619038  39.732599
## 
## $Quantile
##        0%        5%       25%       50%       75%       95%      100% 
##  27.61606  62.61550 107.56766 139.51216 155.21306 189.49942 206.94153 
## 
## $`Confidence interval`
## [1] 120.5449 142.5715
mean_density(wtp_rrm, wtp_rum)
## [1] 0.7308
pnorm( (mean(wtp_rrm)-mean(wtp_rum)) / sqrt(var(wtp_rrm)+var(wtp_rum)), lower.tail = FALSE)
## [1] 0.2914486

Normality test

shapiro.test(wtp_rum)
## 
##  Shapiro-Wilk normality test
## 
## data:  wtp_rum
## W = 0.9558, p-value = 0.0596
shapiro.test(wtp_rrm)
## 
##  Shapiro-Wilk normality test
## 
## data:  wtp_rrm
## W = 0.9641, p-value = 0.1324

welch’s t.test

t_test_general = t.test(wtp_rrm, wtp_rum, var.equal=FALSE)
print(t_test_general)
## 
##  Welch Two Sample t-test
## 
## data:  wtp_rrm and wtp_rum
## t = 3.8831, df = 91.807, p-value = 0.0001942
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  13.43082 41.55488
## sample estimates:
## mean of x mean of y 
##  131.5582  104.0654
cond <- ifelse(t_test_general$p.value <= 0.05, 
               'Reject H0: There is empirical evidence of significant differences in Means of RRM and RUM WTP', 'Accept H0: The means are equal')
cat(sprintf('Test Result: %s', cond))
## Test Result: Reject H0: There is empirical evidence of significant differences in Means of RRM and RUM WTP