RUM - RRM Gender results
marbel
26 de mayo de 2015
Load the results of bootstrapping RUM and RRM models and prepare the data
rm(list = ls(all = TRUE))
lkp_income <- read.table(text='Household_Income
"Less than $25,000"
"$25,000 - $50,000"
"$51,000 - $75,000"
"More than $75,000"',
header=TRUE)library(ggplot2)
library(reshape2)
### RUM general results
load('RUM_wtp_income_opt_out.RData')
variable_levels <- 1:4
rum_wtp <- sapply(variable_levels, function(lv){
sapply(1:50, function(x) RUM_wtp_income_opt_out[[lv]][[x]][['wtp']])
})
rum_wtp <- data.frame(rum_wtp)
names(rum_wtp) <- lkp_income[,1]
### RRM Gender Results
load('RRM_wtp_income_opt_out.RData')
length(RRM_wtp_income_opt_out)## [1] 4rrm_wtp <- sapply(variable_levels, function(gdr){
sapply(1:50, function(x) mean(RRM_wtp_income_opt_out[[gdr]][[x]]))
})
rrm_wtp <- data.frame(rrm_wtp)
names(rrm_wtp) <- lkp_income[,1]Data reshaping to long format for plotting
m_rrm_wtp <- melt(rrm_wtp)## No id variables; using all as measure variablesm_rum_wtp <- melt(rum_wtp)## No id variables; using all as measure variablesm_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('Income', 'Mean_WTP', 'Model')
mdata$Model <- as.factor(mdata$Model)
mdata$Income <- as.factor(mdata$Income)p <- ggplot(mdata, aes(x=Model, y=Mean_WTP, fill=Model))
p <- p + geom_boxplot()
p <- p + theme_bw()
p <- p + facet_wrap(~Income, scales='free')
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Income",
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(~Income, scales='free')
p <- p + labs(list(title="Comparison between RRM and RUM model's \naverage willingness to pay by Income",
x = 'Bootstrapped Willingness to pay', y='Density'))
print(suppressWarnings(p))
distribution comparison
var_levels <- unique(as.character(mdata$Income))
wtp_rum_1 <- sort(mdata[mdata$Model == 'RUM' & mdata$Income %in% var_levels[[1]], 'Mean_WTP'])
wtp_rum_2 <- sort(mdata[mdata$Model == 'RUM' & mdata$Income %in% var_levels[[2]], 'Mean_WTP'])
wtp_rum_3 <- sort(mdata[mdata$Model == 'RUM' & mdata$Income == var_levels[[3]], 'Mean_WTP'])
wtp_rum_4 <- sort(mdata[mdata$Model == 'RUM' & mdata$Income == var_levels[[4]], 'Mean_WTP'])
wtp_rrm_1 <- sort(mdata[mdata$Model == 'RRM' & mdata$Income == var_levels[[1]], 'Mean_WTP'])
wtp_rrm_2 <- sort(mdata[mdata$Model == 'RRM' & mdata$Income == var_levels[[2]], 'Mean_WTP'])
wtp_rrm_3 <- sort(mdata[mdata$Model == 'RRM' & mdata$Income == var_levels[[3]], 'Mean_WTP'])
wtp_rrm_4 <- sort(mdata[mdata$Model == 'RRM' & mdata$Income == var_levels[[4]], 'Mean_WTP'])KS and Wilcox tests comparing each Income data
Income 1
ks.test(wtp_rrm_1, wtp_rum_1)##
## Two-sample Kolmogorov-Smirnov test
##
## data: wtp_rrm_1 and wtp_rum_1
## D = 0.56, p-value = 1.453e-07
## alternative hypothesis: two-sidedwilcox.test(wtp_rrm_1, wtp_rum_1)##
## Wilcoxon rank sum test with continuity correction
##
## data: wtp_rrm_1 and wtp_rum_1
## W = 2082, p-value = 9.913e-09
## alternative hypothesis: true location shift is not equal to 0Income 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-sidedwilcox.test(wtp_rrm_2, wtp_rum_2)##
## Wilcoxon rank sum test with continuity correction
##
## data: wtp_rrm_2 and wtp_rum_2
## W = 2453, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Income 3
ks.test(wtp_rrm_3, wtp_rum_3)##
## Two-sample Kolmogorov-Smirnov test
##
## data: wtp_rrm_3 and wtp_rum_3
## D = 0.72, p-value = 8.756e-13
## alternative hypothesis: two-sidedwilcox.test(wtp_rrm_3, wtp_rum_3)##
## Wilcoxon rank sum test with continuity correction
##
## data: wtp_rrm_3 and wtp_rum_3
## W = 2234, p-value = 1.201e-11
## alternative hypothesis: true location shift is not equal to 0Income 4
ks.test(wtp_rrm_4, wtp_rum_4)##
## Two-sample Kolmogorov-Smirnov test
##
## data: wtp_rrm_4 and wtp_rum_4
## D = 0.82, p-value < 2.2e-16
## alternative hypothesis: two-sidedwilcox.test(wtp_rrm_4, wtp_rum_4)##
## Wilcoxon rank sum test with continuity correction
##
## data: wtp_rrm_4 and wtp_rum_4
## W = 2417, p-value = 8.864e-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] 112.43864 63.98176
##
## $Quantile
## 0% 25% 50% 75% 100%
## -85.42003 86.99779 113.04625 153.75385 244.12507sumarizar(wtp_rum_1)## $`Mean and SD`
## [1] 32.05185 57.20353
##
## $Quantile
## 0% 25% 50% 75% 100%
## -116.472373 -1.883046 30.169910 72.933550 156.849098sumarizar(wtp_rrm_2)## $`Mean and SD`
## [1] 281.05069 68.79175
##
## $Quantile
## 0% 25% 50% 75% 100%
## 142.1670 223.1306 298.1621 333.8635 420.7062sumarizar(wtp_rum_2)## $`Mean and SD`
## [1] 107.61437 53.96365
##
## $Quantile
## 0% 25% 50% 75% 100%
## 0.1102543 68.0576091 114.9252077 147.9097676 189.1568639sumarizar(wtp_rrm_3)## $`Mean and SD`
## [1] 435.2968 153.1868
##
## $Quantile
## 0% 25% 50% 75% 100%
## 176.7179 335.9343 410.2086 539.9644 885.3688sumarizar(wtp_rum_3)## $`Mean and SD`
## [1] 225.85238 91.63732
##
## $Quantile
## 0% 25% 50% 75% 100%
## 12.00538 172.24621 231.36733 282.07493 534.09528sumarizar(wtp_rrm_4)## $`Mean and SD`
## [1] 285.67524 94.98018
##
## $Quantile
## 0% 25% 50% 75% 100%
## 94.89067 231.33590 279.10056 338.83628 561.60220sumarizar(wtp_rum_4)## $`Mean and SD`
## [1] 103.40865 49.51403
##
## $Quantile
## 0% 25% 50% 75% 100%
## -13.28001 75.66568 104.09995 130.35390 229.07225mean_density(wtp_rum_1, wtp_rrm_1)## [1] 0.1672mean_density(wtp_rum_2, wtp_rrm_2)## [1] 0.0188mean_density(wtp_rum_3, wtp_rrm_3)## [1] 0.1064mean_density(wtp_rum_4, wtp_rrm_4)## [1] 0.0332pnorm( (mean(wtp_rrm_1)-mean(wtp_rum_1)) /
sqrt(var(wtp_rrm_1) + var(wtp_rum_1)), lower.tail = FALSE)## [1] 0.1744725pnorm( (mean(wtp_rrm_2)-mean(wtp_rum_2)) /
sqrt(var(wtp_rrm_2) + var(wtp_rum_2)), lower.tail = FALSE)## [1] 0.02364645pnorm( (mean(wtp_rrm_3)-mean(wtp_rum_3)) /
sqrt(var(wtp_rrm_3) + var(wtp_rum_3)), lower.tail = FALSE)## [1] 0.1203311pnorm( (mean(wtp_rrm_4)-mean(wtp_rum_4)) /
sqrt(var(wtp_rrm_4) + var(wtp_rum_4)), lower.tail = FALSE)## [1] 0.04441022