rm(list = ls(all = TRUE))
lkp <- c('Male', 'Female')

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

library(ggplot2)
library(reshape2)
### RUM general results
load('RUM_wtp_gender_opt_out.RData')

rum_wtp <- sapply(1:2, function(gdr){
  abs(sapply(1:50, function(x) RUM_wtp_gender_opt_out[[gdr]][[x]][['wtp']]))
})
rum_wtp <- data.frame(rum_wtp)
names(rum_wtp) <- lkp

### RRM Gender Results
load('RRM_wtp_gender_opt_out.RData')
length(RRM_wtp_gender_opt_out)
## [1] 2
rrm_wtp <- sapply(1:2, function(gdr){
  abs(sapply(1:50, function(x) mean(RRM_wtp_gender_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('Gender', 'Mean_WTP', 'Model')
mdata$Model <- as.factor(mdata$Model)
mdata$Gender <- as.factor(mdata$Gender)
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 + facet_wrap(~Gender, scales='free')
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Gender",
                   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 + theme_bw()
p <- p + scale_fill_manual(values=c('grey40', 'grey70'))
p <- p + facet_wrap(~Gender, scales='free')
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Gender",
                   x = 'Bootstrapped Willingness to pay', y='Density'))
print(suppressWarnings(p))

distribution comparison

gender_1 = 'Male'
gender_2 = 'Female'

wtp_rum_gender_1 <- sort(mdata[mdata$Model == 'RUM' & mdata$Gender == gender_1, 'Mean_WTP'])
wtp_rum_gender_2 <- sort(mdata[mdata$Model == 'RUM' & mdata$Gender == gender_2, 'Mean_WTP'])

wtp_rrm_gender_1 <- sort(mdata[mdata$Model == 'RRM' & mdata$Gender == gender_1, 'Mean_WTP'])
wtp_rrm_gender_2 <- sort(mdata[mdata$Model == 'RRM' & mdata$Gender == gender_2, 'Mean_WTP'])

KS and Wilcox tests comparing each gender data

Gender 1 - Male

ks.test(wtp_rrm_gender_1, wtp_rum_gender_1)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_gender_1 and wtp_rum_gender_1
## D = 0.38, p-value = 0.001315
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm_gender_1, wtp_rum_gender_1)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm_gender_1 and wtp_rum_gender_1
## W = 1785, p-value = 0.0002289
## alternative hypothesis: true location shift is not equal to 0

Gender 2 - Female

ks.test(wtp_rrm_gender_2, wtp_rum_gender_2)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  wtp_rrm_gender_2 and wtp_rum_gender_2
## D = 0.24, p-value = 0.1124
## alternative hypothesis: two-sided
wilcox.test(wtp_rrm_gender_2, wtp_rum_gender_2)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  wtp_rrm_gender_2 and wtp_rum_gender_2
## W = 1521, p-value = 0.06221
## alternative hypothesis: true location shift is not equal to 0
library('plotrix')
## Warning: package 'plotrix' was built under R version 3.1.3
# 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'
  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_gender_1)
## $`Mean - SE - SD`
## [1] 185.858009   5.547762  39.228601
## 
## $Quantile
##        0%        5%       25%       50%       75%       95%      100% 
##  86.43856 118.60536 160.07508 191.32628 203.72643 238.83334 283.10770
sumarizar(wtp_rum_gender_1)
## $`Mean - SE - SD`
## [1] 156.866245   5.226988  36.960387
## 
## $Quantile
##        0%        5%       25%       50%       75%       95%      100% 
##  68.96867  95.71418 133.15079 158.86953 176.53708 208.84475 258.33554
sumarizar(wtp_rrm_gender_2)
## $`Mean - SE - SD`
## [1] 91.374395  6.404918 45.289608
## 
## $Quantile
##         0%         5%        25%        50%        75%        95% 
##   5.933242  23.402361  64.500333  84.486250 124.977258 157.323386 
##       100% 
## 234.473014
sumarizar(wtp_rum_gender_2)
## $`Mean - SE - SD`
## [1] 76.007783  5.409549 38.251288
## 
## $Quantile
##         0%         5%        25%        50%        75%        95% 
##   4.912536  18.288110  53.238597  70.768017 103.449154 130.415845 
##       100% 
## 192.612729
mean_density(wtp_rrm_gender_1, wtp_rum_gender_1)
## [1] 0.714
mean_density(wtp_rrm_gender_2, wtp_rum_gender_2)
## [1] 0.6084
pnorm( (mean(wtp_rrm_gender_1)-mean(wtp_rum_gender_1)) / 
         sqrt(var(wtp_rrm_gender_1) + var(wtp_rum_gender_1)), lower.tail = FALSE)
## [1] 0.2953215
pnorm( (mean(wtp_rrm_gender_2)-mean(wtp_rum_gender_2)) / 
         sqrt(var(wtp_rrm_gender_2) + var(wtp_rum_gender_2)), lower.tail = FALSE)
## [1] 0.3977351

Normality test

shapiro.test(wtp_rrm_gender_1)
## 
##  Shapiro-Wilk normality test
## 
## data:  wtp_rrm_gender_1
## W = 0.9857, p-value = 0.8029
shapiro.test(wtp_rum_gender_1)
## 
##  Shapiro-Wilk normality test
## 
## data:  wtp_rum_gender_1
## W = 0.9891, p-value = 0.922
shapiro.test(wtp_rrm_gender_2)
## 
##  Shapiro-Wilk normality test
## 
## data:  wtp_rrm_gender_2
## W = 0.9689, p-value = 0.2092
shapiro.test(wtp_rum_gender_2)
## 
##  Shapiro-Wilk normality test
## 
## data:  wtp_rum_gender_2
## W = 0.9689, p-value = 0.2093

One sided t.test

t.test(wtp_rrm_gender_1 - wtp_rum_gender_1, mu=0, var.equal=FALSE, alternative='greater')
## 
##  One Sample t-test
## 
## data:  wtp_rrm_gender_1 - wtp_rum_gender_1
## t = 54.9369, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  28.107    Inf
## sample estimates:
## mean of x 
##  28.99176
boxplot(wtp_rrm_gender_1 - wtp_rum_gender_1)

t.test(wtp_rrm_gender_2 - wtp_rum_gender_2, mu=0, var.equal=FALSE, alternative='greater')
## 
##  One Sample t-test
## 
## data:  wtp_rrm_gender_2 - wtp_rum_gender_2
## t = 15.0045, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  13.6496     Inf
## sample estimates:
## mean of x 
##  15.36661
boxplot(wtp_rrm_gender_2 - wtp_rum_gender_2)

Welch’s Two Sample t-test

Alternative Greater

t_test_gender_1 = t.test(wtp_rrm_gender_1, wtp_rum_gender_1, 
                         var.equal=FALSE)
print(t_test_gender_1)
## 
##  Welch Two Sample t-test
## 
## data:  wtp_rrm_gender_1 and wtp_rum_gender_1
## t = 3.8036, df = 97.654, p-value = 0.0002486
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  13.86494 44.11859
## sample estimates:
## mean of x mean of y 
##  185.8580  156.8662
t_test_gender_2 = t.test(wtp_rrm_gender_2, wtp_rum_gender_2, 
                         var.equal=FALSE)
print(t_test_gender_2)
## 
##  Welch Two Sample t-test
## 
## data:  wtp_rrm_gender_2 and wtp_rum_gender_2
## t = 1.8329, df = 95.331, p-value = 0.06994
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.276363 32.009587
## sample estimates:
## mean of x mean of y 
##  91.37439  76.00778
# H0 Means difference of WTP_RUM - WTP_RRM are equal to 0
# H1 Means difference WTP_RUM - WTP_RRM are not equal to 0

test_summary <- function(t_test) {
  res = ifelse(t_test$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
p-value: %s', res, t_test$p.value))

}

welch’s test results

test_summary(t_test_gender_1)
## Test Result: 
## Reject H0: There is empirical evidence of significant differences in Means of RRM and RUM WTP
## p-value: 0.00024860332192103
test_summary(t_test_gender_2)
## Test Result: 
## Accept H0: The means are equal
## p-value: 0.0699358577788875