Load the results of bootstrapping RUM and RRM models and prepare the data

rm(list = ls(all = TRUE))
lkp <- c('18-35', '36-55', '> 55')
library(ggplot2)
library(reshape2)
### RUM Age results
load('RUM_wtp_age_opt_out.RData')

variable_levels <- 1:3
rum_wtp <- sapply(variable_levels, function(lv){
  sapply(1:50, function(x) RUM_wtp_age_opt_out[[lv]][[x]][['wtp']])
})
rum_wtp <- data.frame(rum_wtp)
names(rum_wtp) <- lkp

### RRM Age Results
load('RRM_wtp_age_opt_out.RData')
length(RRM_wtp_age_opt_out)
## [1] 3
rrm_wtp <- sapply(variable_levels, function(gdr){
  sapply(1:50, function(x) mean(RRM_wtp_age_opt_out[[gdr]][[x]]))
})
rrm_wtp <- data.frame(rrm_wtp)
names(rrm_wtp) <- lkp

Data reshaping to long format for plotting

m_rrm_wtp <- melt(rrm_wtp)
## No id variables; using all as measure variables
m_rum_wtp <- melt(rum_wtp)
## No id variables; using all as measure variables
m_rrm_wtp$Model <- 'RRM'
m_rum_wtp$Model <- 'RUM'

# Prepare the data for the plot
mdata <- data.frame(rbind(m_rrm_wtp, m_rum_wtp))
names(mdata) <- c('Age', 'Mean_WTP', 'Model')
mdata$Model <- as.factor(mdata$Model)
mdata$Age <- as.factor(mdata$Age)
p <- ggplot(mdata, aes(x=Model, y=Mean_WTP, fill=Model))
p <- p + geom_boxplot()
p <- p + theme_bw()
p <- p + facet_wrap(~Age)
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Age",
                   y = 'Average willingness to pay'))
print(suppressWarnings(p))

p <- ggplot(mdata, aes(x=Mean_WTP, fill=Model))
p <- p + geom_histogram(aes(y = ..density..), alpha=.5) + geom_density(alpha=0.2)
p <- p + theme_bw()
p <- p + facet_wrap(~Age)
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Age",
                   x = 'Bootstrapped Willingness to pay', y='Density'))
print(suppressWarnings(p))

Distribution comparison

var_levels <- unique(as.character(mdata$Age))

wtp_rum_1 <- sort(mdata[mdata$Model == 'RUM' & mdata$Age %in% var_levels[[1]], 'Mean_WTP'])
wtp_rum_2 <- sort(mdata[mdata$Model == 'RUM' & mdata$Age %in% var_levels[[2]], 'Mean_WTP'])
wtp_rum_3 <- sort(mdata[mdata$Model == 'RUM' & mdata$Age == var_levels[[3]], 'Mean_WTP'])

wtp_rrm_1 <- sort(mdata[mdata$Model == 'RRM' & mdata$Age == var_levels[[1]], 'Mean_WTP'])
wtp_rrm_2 <- sort(mdata[mdata$Model == 'RRM' & mdata$Age == var_levels[[2]], 'Mean_WTP'])
wtp_rrm_3 <- sort(mdata[mdata$Model == 'RRM' & mdata$Age == var_levels[[3]], 'Mean_WTP'])

KS and Wilcox tests comparing each Age data

Age 1

ks.test(wtp_rrm_1, wtp_rum_1)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_1 and wtp_rum_1
## D = 0.58, p-value = 4.048e-08
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm_1, wtp_rum_1)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm_1 and wtp_rum_1
## W = 2089, p-value = 7.449e-09
## alternative hypothesis: true location shift is not equal to 0

Age 2

ks.test(wtp_rrm_2, wtp_rum_2)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_2 and wtp_rum_2
## D = 0.9, p-value < 2.2e-16
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm_2, wtp_rum_2)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm_2 and wtp_rum_2
## W = 2423, p-value = 6.319e-16
## alternative hypothesis: true location shift is not equal to 0

Age 3

ks.test(wtp_rrm_3, wtp_rum_3)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_3 and wtp_rum_3
## D = 0.94, p-value < 2.2e-16
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm_3, wtp_rum_3)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm_3 and wtp_rum_3
## W = 2483, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
# summary functions used below
sumarizar <- function(x){
  ls = list(c(mean(x), sd(x)), quantile(x))
  names(ls)[[1]] <- 'Mean and SD'
  names(ls)[[2]] <- 'Quantile'
  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_rrm_1)
## $`Mean and SD`
## [1] 136.924  72.748
## 
## $Quantile
##        0%       25%       50%       75%      100% 
## -41.10289  95.52695 144.03033 181.59968 300.33235
sumarizar(wtp_rum_1)
## $`Mean and SD`
## [1] 48.25875 60.46688
## 
## $Quantile
##           0%          25%          50%          75%         100% 
## -104.8547865    0.5578088   58.2334902   83.6093036  185.1472378
sumarizar(wtp_rrm_2)
## $`Mean and SD`
## [1] 226.73150  57.99228
## 
## $Quantile
##        0%       25%       50%       75%      100% 
##  75.85496 180.54285 222.84887 273.68463 330.31393
sumarizar(wtp_rum_2)
## $`Mean and SD`
## [1] 95.94442 41.21492
## 
## $Quantile
##         0%        25%        50%        75%       100% 
##  -6.542612  66.679660 100.584700 121.005547 170.085405
sumarizar(wtp_rrm_3)
## $`Mean and SD`
## [1] 339.95845  65.13545
## 
## $Quantile
##       0%      25%      50%      75%     100% 
## 196.5012 303.6595 329.3239 367.7186 486.0870
sumarizar(wtp_rum_3)
## $`Mean and SD`
## [1] 149.9936  45.8120
## 
## $Quantile
##        0%       25%       50%       75%      100% 
##  37.22631 126.95968 143.33617 167.80172 245.82271
mean_density(wtp_rum_1, wtp_rrm_1)
## [1] 0.1644
mean_density(wtp_rum_2, wtp_rrm_2)
## [1] 0.0308
mean_density(wtp_rum_3, wtp_rrm_3)
## [1] 0.0068
pnorm( (mean(wtp_rrm_1)-mean(wtp_rum_1)) / 
         sqrt(var(wtp_rrm_1) + var(wtp_rum_1)), lower.tail = FALSE)
## [1] 0.1743025
pnorm( (mean(wtp_rrm_2)-mean(wtp_rum_2)) / 
         sqrt(var(wtp_rrm_2) + var(wtp_rum_2)), lower.tail = FALSE)
## [1] 0.03301003
pnorm( (mean(wtp_rrm_3)-mean(wtp_rum_3)) / 
         sqrt(var(wtp_rrm_3) + var(wtp_rum_3)), lower.tail = FALSE)
## [1] 0.008527652