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

rm(list = ls(all = TRUE))
lkp <- c('White', 'Black', 'Indian, Alaskan Native, Asian, Other', 'Hispanic or Latino')
print(cbind(1:4, lkp))
##          lkp                                   
## [1,] "1" "White"                               
## [2,] "2" "Black"                               
## [3,] "3" "Indian, Alaskan Native, Asian, Other"
## [4,] "4" "Hispanic or Latino"
library(ggplot2)
library(reshape2)
### RUM general results
load('RUM_wtp_ethnic_opt_out.RData')

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

### RRM Gender Results
load('RRM_wtp_ethnic_opt_out.RData')
length(RRM_wtp_ethnic_opt_out)
## [1] 4
rrm_wtp <- sapply(variable_levels, function(gdr){
  sapply(1:50, function(x) mean(RRM_wtp_ethnic_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('Ethnicity', 'Mean_WTP', 'Model')
mdata$Model <- as.factor(mdata$Model)
mdata$Ethnicity <- as.factor(mdata$Ethnicity)
p <- ggplot(mdata, aes(x=Model, y=Mean_WTP, fill=Model))
p <- p + geom_boxplot()
p <- p + theme_bw()
p <- p + facet_wrap(~Ethnicity, scales='free')
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Ethnicity",
                   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(~Ethnicity, scales='free')
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Ethnicity",
                   x = 'Bootstrapped Willingness to pay', y='Density'))
print(suppressWarnings(p))

distribution comparison

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

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

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

KS and Wilcox tests comparing each Income data

Ethnicity 1 - White

ks.test(wtp_rrm_1, wtp_rum_1)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_1 and wtp_rum_1
## D = 0.92, p-value < 2.2e-16
## 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 = 2470, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Ethnicity 2 - Black

ks.test(wtp_rrm_2, wtp_rum_2)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_2 and wtp_rum_2
## D = 0.56, p-value = 1.453e-07
## 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 = 2074, p-value = 1.37e-08
## alternative hypothesis: true location shift is not equal to 0

Ethnicity 3 - Indian or Alaskan Native, Asian, Other

ks.test(wtp_rrm_3, wtp_rum_3)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_3 and wtp_rum_3
## D = 0.88, 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 = 2358, p-value = 2.26e-14
## alternative hypothesis: true location shift is not equal to 0

Ethnicity 4 - Hispanic or Latino

ks.test(wtp_rrm_4, wtp_rum_4)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_4 and wtp_rum_4
## D = 0.4, p-value = 0.0005823
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm_4, wtp_rum_4)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm_4 and wtp_rum_4
## W = 1749, p-value = 0.0005891
## 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)
}

pnorm2 <- function(wtp_rum_1, wtp_rrm_1){
  pnorm( (mean(wtp_rrm_1)-mean(wtp_rum_1)) / 
         sqrt(var(wtp_rrm_1) + var(wtp_rum_1)), lower.tail = FALSE)
}

Rest of the summary results

sumarizar(wtp_rrm_1)
## $`Mean and SD`
## [1] 210.8227  47.8646
## 
## $Quantile
##        0%       25%       50%       75%      100% 
##  97.75013 176.65387 210.38585 241.38409 318.66766
sumarizar(wtp_rum_1)
## $`Mean and SD`
## [1] 87.10588 31.82039
## 
## $Quantile
##        0%       25%       50%       75%      100% 
##  11.64370  59.99135  85.45423 107.59886 157.26874
sumarizar(wtp_rrm_2)
## $`Mean and SD`
## [1] 309.3654 150.0753
## 
## $Quantile
##        0%       25%       50%       75%      100% 
## -108.8755  227.5405  334.2699  397.4069  634.0124
sumarizar(wtp_rum_2)
## $`Mean and SD`
## [1] 130.2062 114.6462
## 
## $Quantile
##         0%        25%        50%        75%       100% 
## -182.27912   72.97888  122.32012  213.80035  334.69260
sumarizar(wtp_rrm_3)
## $`Mean and SD`
## [1] 1238.030 1322.068
## 
## $Quantile
##         0%        25%        50%        75%       100% 
##   66.33406  501.77822  807.83385 1277.19687 7088.83229
sumarizar(wtp_rum_3)
## $`Mean and SD`
## [1] 200.8518 140.4476
## 
## $Quantile
##       0%      25%      50%      75%     100% 
## -69.9917  89.7675 210.3698 306.3625 547.4369
sumarizar(wtp_rrm_4)
## $`Mean and SD`
## [1] 294.5571 111.6429
## 
## $Quantile
##       0%      25%      50%      75%     100% 
##  40.6505 232.1739 297.3461 382.3721 539.2188
sumarizar(wtp_rum_4)
## $`Mean and SD`
## [1] 210.4596 126.1111
## 
## $Quantile
##           0%          25%          50%          75%         100% 
##   0.05608765 143.81769296 221.66841932 282.46722414 492.63704548
# Mean density
mean_density(wtp_rum_1, wtp_rrm_1)
## [1] 0.012
mean_density(wtp_rum_2, wtp_rrm_2)
## [1] 0.1704
mean_density(wtp_rum_3, wtp_rrm_3)
## [1] 0.0568
mean_density(wtp_rum_4, wtp_rrm_4)
## [1] 0.3004
# pnorm2 is Defined above
pnorm2(wtp_rum_1, wtp_rrm_1)
## [1] 0.01568004
pnorm2(wtp_rum_2, wtp_rrm_2)
## [1] 0.1713972
pnorm2(wtp_rum_3, wtp_rrm_3)
## [1] 0.2176595
pnorm2(wtp_rum_4, wtp_rrm_4)
## [1] 0.3087815