########################################################### ### R code for examples in Chapter 3 (ver.140808) ### ### ### ### Hideo Aizaki, Tomoaki Nakatani, Kazuo Sato (2014) ### ### Stated Preference Methods Using R (SPMUR) ### ### Chapman & Hall/CRC Press ### ### ### ### This file is licensed under the GNU GPLv2 or later ### ### and is available from ### ### http://www.agr.hokudai.ac.jp/spmur/ ### ### ### ### NOTE: Some lines of code in this file are not shown ### ### in the book explicitly. A comment "not shown ### ### in SPMUR" is added to such lines. ### ########################################################### ########################################################### ### 3.3.1 Overview ########################################################### library(support.CEs) library(survival) library(mded) ########################################################### ### 3.3.2 Creating a DCE design ########################################################### rd <- rotation.design( attribute.names = list(X = c("x1", "x2", "x3"), Y = c("y1", "y2", "y3"), Z = c("10", "20", "30")), nalternatives = 2, nblocks = 3, seed = 123) rd rd$alternatives$alt.1 ########################################################### ### 3.3.3 Converting a DCE design into questionnaire format ########################################################### questionnaire(choice.experiment.design = rd) ########################################################### ### 3.3.4 Creating a design matrix ########################################################### dm.rd <- make.design.matrix( choice.experiment.design = rd, optout = TRUE, categorical.attributes = c("X", "Y"), continuous.attributes = c("Z"), unlabeled = TRUE, common = NULL, binary = FALSE) dm.rd ########################################################### ### 3.3.5 Creating a data set ########################################################### res <- data.frame(ID = c(1:30), BLOCK = rep(c(1:3), times = 10), q1 = c(2, 2, 2, 1, 2, 1, 1, 2, 1, 3, 2, 1, 3, 2, 1, 2, 2, 1, 1, 2, 1, 1, 2, 1, 3, 2, 1, 3, 2, 1), q2 = c(1, 2, 3, 1, 2, 2, 1, 2, 1, 3, 2, 2, 3, 2, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 1, 2, 1, 3, 2, 1), q3 = c(2, 1, 2, 1, 1, 2, 1, 3, 1, 1, 1, 2, 1, 3, 2, 1, 1, 2, 3, 2, 3, 3, 2, 1, 1, 2, 2, 1, 3, 3)) ds <- make.dataset(respondent.dataset = res, design.matrix = dm.rd, choice.indicators = c("q1", "q2", "q3")) ds[1:10, ] ds[268:270, ] ########################################################### ### 3.3.6 Conducting statistical analysis ########################################################### fm <- RES ~ ASC + x2 + x3 + y2 + y3 + Z + strata(STR) cl <- clogit(fm, data = ds) cl cl$coefficients cl$var ########################################################### ### 3.3.7 Calculating goodness-of-fit measures ########################################################### gofm(cl) ########################################################### ### 3.3.8 Calculating MWTPs ########################################################### mwtp.cl <- mwtp(output = cl, monetary.variables=c("Z"), nrep = 1000, confidence.level = 0.95, method = "kr", seed = 123) mwtp.cl ########################################################### ### 3.3.9 Testing the difference between two independent ### MWTP distributions ########################################################### MWTP_x2 <- mwtp.cl$mwtps[, "x2"] MWTP_x3 <- mwtp.cl$mwtps[, "x3"] mded.cl <- mded(distr1 = MWTP_x3, distr2 = MWTP_x2) mded.cl ########################################################### ### 3.4.1 Unlabeled DCE example ########################################################### d.rice <- rotation.design( attribute.names = list( Region = c("RegA", "RegB", "RegC"), Cultivation = c("Conv", "NoChem", "Organic"), Price = c("1700", "2000", "2300")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, # mix-and-match method seed = 987) d.rice questionnaire(choice.experiment.design = d.rice) data(rice) head(rice) dm.rice <- make.design.matrix( choice.experiment.design = d.rice, optout = TRUE, # include opt-out option categorical.attributes = c("Region", "Cultivation"), continuous.attributes = c("Price"), unlabeled = TRUE) # unlabeled design head(dm.rice) ds.rice <- make.dataset( respondent.dataset = rice, choice.indicators = c("q1", "q2", "q3", "q4", "q5", "q6", "q7", "q8", "q9"), design.matrix = dm.rice) head(ds.rice) fm.rice <- RES ~ ASC + RegB + RegC + NoChem + Organic + NoChem:F + Organic:F + Price + strata(STR) out.rice <- clogit(fm.rice, data = ds.rice) out.rice gofm(out.rice) out.rice.r <- update(out.rice, .~. - RegB - RegC + I(RegB + RegC)) -2 * (out.rice.r$loglik[2] - out.rice$loglik[2]) mwtp(output = out.rice, monetary.variables = c("Price"), nonmonetary.variables = c("RegB", "RegC", "NoChem", "Organic", "NoChem:F", "Organic:F"), confidence.level = 0.95, method = "kr", seed = 987) ########################################################### ### 3.4.2 Labeled design example ########################################################### d.pork <- Lma.design( attribute.names = list( Price = c("100", "130", "160", "190")), nalternatives = 3, nblocks = 4, row.renames = FALSE, seed = 987) d.pork questionnaire(choice.experiment.design = d.pork) data(pork) head(pork) dm.pork <- make.design.matrix( choice.experiment.design = d.pork, optout = TRUE, # include opt-out option continuous.attributes = c("Price"), unlabeled = FALSE) # labeled design ds.pork <- make.dataset( respondent.dataset = pork, choice.indicators = c("q1", "q2", "q3", "q4"), design.matrix = dm.pork) fm.pork <- RES ~ ASC1 + Price1 + # variables for Imported ASC2 + Price2 + # variables for Domestic ASC3 + Price3 + # variables for HQ domestic strata(STR) # stratification variable out.pork <- clogit(fm.pork, data = ds.pork) out.pork gofm(out.pork) mwtp.pork <- mwtp( output = out.pork, monetary.variables = c("Price1", "Price2", "Price3"), nonmonetary.variables = list( c("ASC1"), c("ASC2"), c("ASC3")), confidence.level = 0.95, method = "kr", seed = 987) mwtp.pork hist(mwtp.pork$mwtps[, "ASC1"], main = "", xlab = "ASC1") hist(mwtp.pork$mwtps[, "ASC2"], main = "", xlab = "ASC2") hist(mwtp.pork$mwtps[, "ASC3"], main = "", xlab = "ASC3") ds.pork$Price <- ds.pork$Price1 + ds.pork$Price2 + ds.pork$Price3 ds.pork[1:4, c("Price1", "Price2", "Price3", "Price")] fm.porkg <- RES ~ ASC1 + ASC2 + ASC3 + # ASC for each pork Price + # generic price variable strata(STR) # stratification variable out.porkg <- clogit(fm.porkg, data = ds.pork) out.porkg gofm(out.porkg) out.pork$loglik out.porkg$loglik -2 * (out.porkg$loglik[2] - out.pork$loglik[2]) Price.imp <- c(135) Price.dom <- c(150) Price.hqd <- c(100:190) out.pork$coef v.imp <- out.pork$coef[1] + out.pork$coef[2] * Price.imp v.dom <- out.pork$coef[3] + out.pork$coef[4] * Price.dom v.hqd <- out.pork$coef[5] + out.pork$coef[6] * Price.hqd denominator.pork <- exp(v.imp) + exp(v.dom) + exp(v.hqd) + 1 Prob.imp <- exp(v.imp) / denominator.pork Prob.dom <- exp(v.dom) / denominator.pork Prob.hqd <- exp(v.hqd) / denominator.pork Prob.not <- 1 - (Prob.imp + Prob.dom + Prob.hqd) plot(x = Price.hqd, # x-axis y = Prob.imp, # y-axis ylab = "Choice probability", # label for y-axis ylim = c(0, 1), # limits for y-axis type = "l", # plot type is line lty = "solid") # line type is solid line lines(Price.hqd, Prob.dom, # add dashed line of Prob.dom lty = "dashed") # on plot lines(Price.hqd, Prob.hqd, # add dotted line of Prob.hqd lty = "dotted") # on plot lines(Price.hqd, Prob.not, # add two-dashed line of Prob.not lty = "twodash") # on plot legend("topright", # add legend at top right of plot legend = c("Imported pork", # set names of each item "Domestic pork", "HQ domestic pork", "None-of-these"), lty = c("solid", # set line types of each item "dashed", # in the order corresponding to "dotted", # that of the item names "twodash")) ########################################################### ### 3.4.3 BDCE example ########################################################### d.rural <- Lma.design( attribute.names = list( Area = c("20", "40", "60", "80"), Facility = c("None", "Agr", "Env", "Rec"), Tax = c("1000", "3000", "5000", "7000")), nalternatives = 1, # create design for a BDCE nblocks = 4, row.renames = FALSE, seed = 987) d.rural common.alt <- c(Area = "0", Facility = "None", Tax = "0") questionnaire(choice.experiment.design = d.rural, common = common.alt) data(rural) # extract rows with Region = 1 rural1 <- subset(rural, Region == 1) # extract rows with Region = 2 rural2 <- subset(rural, Region == 2) head(rural1) head(rural2) dm.rural <- make.design.matrix( choice.experiment.design = d.rural, optout = FALSE, # do not include opt-out option categorical.attributes = c("Facility"), continuous.attributes = c("Area", "Tax"), unlabeled = TRUE, # unlabeled design common = common.alt, # include common alternative binary = TRUE) # BDCEs head(dm.rural) dm.ruralop <- make.design.matrix( choice.experiment.design = d.rural, optout = TRUE, # include opt-out option categorical.attributes = c("Facility"), continuous.attributes = c("Area", "Tax"), unlabeled = TRUE, # unlabeled design common = NULL, # do not include common alternative binary = TRUE) # BDCEs # test whether dm.rural is equal to dm.ruralop identical(dm.rural, dm.ruralop) ds.rural1 <- make.dataset( respondent.dataset = rural1, design.matrix = dm.rural, choice.indicators = c("q1", "q2", "q3", "q4"), detail = FALSE) ds.rural2 <- make.dataset( respondent.dataset = rural2, design.matrix = dm.rural, choice.indicators = c("q1", "q2", "q3", "q4"), detail = FALSE) head(ds.rural1) fm.rural <- RES ~ Agr + Env + Rec + Area + Tax out.rural1 <- glm(fm.rural, family = binomial(link = "logit"), data = ds.rural1) summary(out.rural1) gofm(out.rural1) out.rural2 <- glm(fm.rural, family = binomial(link = "logit"), data = ds.rural2) summary(out.rural2) gofm(out.rural2) mwtp.rural1 <- mwtp(output = out.rural1, monetary.variables = c("Tax"), nonmonetary.variables = c("Agr", "Env", "Rec", "Area"), nreplications = 1000, confidence.level = 0.95, method = "kr", seed = 987) mwtp.rural1 mwtp.rural2 <- mwtp(output = out.rural2, monetary.variables = c("Tax"), nonmonetary.variables = c("Agr", "Env", "Rec", "Area"), nreplications = 1000, confidence.level = 0.95, method = "kr", seed = 987) mwtp.rural2 str(mwtp.rural1$mwtps) MWTP_Rec1 <- mwtp.rural1$mwtps[, "Rec"] MWTP_Rec2 <- mwtp.rural2$mwtps[, "Rec"] mded.rec <- mded(distr1 = MWTP_Rec1, distr2 = MWTP_Rec2, detail = TRUE) mded.rec par(mfrow = c(1, 2)) # not shown in SPMUR hist(mded.rec$diff, main = "(a) Rec", xlab = "Difference") MWTP_Area1 <- mwtp.rural1$mwtps[, "Area"] MWTP_Area2 <- mwtp.rural2$mwtps[, "Area"] mded.area <- mded(distr1 = MWTP_Area1, distr2 = MWTP_Area2, detail = TRUE) mded.area hist(mded.area$diff, main = "(b) Area", xlab = "Difference") v0 <- 0 Plan1 <- c(1, 0, 0, 0, 20, 0) v1 <- sum(out.rural1$coef * Plan1) (-1 / -out.rural1$coef["Tax"]) * (v0 - v1) Plan2 <- c(1, 1, 0, 0, 60, 0) v2 <- sum(out.rural1$coef * Plan2) (-1 / -out.rural1$coef["Tax"]) * (v0 - v2) set.seed(123) COEF <- mvrnorm(out.rural1$coef, # coefficients vcov(out.rural1), # variance-covariance n = 1000) # number of replications V0 <- rep(0, 1000) V1 <- COEF %*% Plan1 # %*% is matrix multiplication C1 <- (-1 / -COEF[, "Tax"]) * (V0 - V1) V2 <- COEF %*% Plan2 C2 <- (-1 / -COEF[, "Tax"]) * (V0 - V2) hist(C1, main = "Case 1", xlab = "Compensating variation") hist(C2, main = "Case 2", xlab = "Compensating variation") quantile(C1, c(0.025, 0.975)) quantile(C2, c(0.025, 0.975)) (-1 / -out.rural1$coef["Tax"]) * (log(exp(v0) + exp(v1)) - log(exp(v0) + exp(v2))) par(mfrow = c(1, 1)) # not shown in SPMUR C3 <- (-1 / -COEF[ ,"Tax"]) * (log(exp(V0) + exp(V1)) - log(exp(V0) + exp(V2))) hist(C3, main = "", xlab = "Compensating variation") quantile(C3, c(0.025, 0.975)) ### END ###