###########################################################
### R code for examples in Chapter 4 (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.             ###
###########################################################


###########################################################
### 4.3.1 Overview
###########################################################
library(DoE.base)
library(crossdes)
library(support.BWS)
library(survival)


###########################################################
### 4.3.2 Constructing choice sets
###########################################################
set.seed(123)
oa <- oa.design(nfactors = 9, nlevels = 2)
oa

set.seed(123)
bibd <- find.BIB(trt = 7, k = 4, b = 7)
bibd
isGYD(bibd)


###########################################################
### 4.3.3 Preparing BWS questions
###########################################################
items <- c("A", "B", "C", "D", "E", "F", "G")
bws.questionnaire(choice.sets = bibd,
                  design.type = 2,  # BIBD
                  item.names = items)


###########################################################
### 4.3.4 Creating the data set
###########################################################
res <- data.frame( # row number format
  ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
  B1 = c(2, 1, 2, 4, 2, 2, 2, 1, 2, 1),
  W1 = c(3, 4, 3, 3, 1, 3, 3, 2, 3, 3),
  B2 = c(4, 3, 3, 3, 3, 1, 1, 2, 1, 1),
  W2 = c(3, 2, 4, 1, 2, 3, 3, 4, 2, 4),
  B3 = c(3, 1, 1, 1, 1, 1, 1, 1, 2, 1),
  W3 = c(1, 4, 2, 2, 4, 4, 2, 3, 3, 3),
  B4 = c(2, 2, 1, 3, 2, 2, 2, 2, 4, 1),
  W4 = c(4, 4, 3, 4, 1, 3, 4, 1, 2, 4),
  B5 = c(1, 3, 2, 1, 3, 2, 1, 1, 1, 1),
  W5 = c(3, 1, 4, 4, 1, 4, 3, 2, 4, 3),
  B6 = c(2, 1, 1, 3, 2, 4, 4, 3, 3, 3),
  W6 = c(3, 2, 3, 4, 3, 2, 2, 4, 4, 2),
  B7 = c(2, 1, 3, 1, 3, 2, 3, 3, 2, 2),
  W7 = c(4, 4, 4, 4, 4, 1, 4, 1, 4, 4))
dat <- bws.dataset(respondent.dataset = res,
                   response.type = 1,  # row number format
                   choice.sets = bibd,
                   design.type = 2,    # BIBD
                   item.names = items)
dat[1:85,]

res2 <- data.frame( # item number format
  ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
  B1 = c(4, 1, 4, 7, 4, 4, 4, 1, 4, 1),
  W1 = c(6, 7, 6, 6, 1, 6, 6, 4, 6, 6),
  B2 = c(6, 4, 4, 4, 4, 2, 2, 3, 2, 2),
  W2 = c(4, 3, 6, 2, 3, 4, 4, 6, 3, 6),
  B3 = c(5, 2, 2, 2, 2, 2, 2, 2, 4, 2),
  W3 = c(2, 7, 4, 4, 7, 7, 4, 5, 5, 5),
  B4 = c(2, 2, 1, 5, 2, 2, 2, 2, 6, 1),
  W4 = c(6, 6, 5, 6, 1, 5, 6, 1, 2, 6),
  B5 = c(1, 4, 3, 1, 4, 3, 1, 1, 1, 1),
  W5 = c(4, 1, 5, 5, 1, 5, 4, 3, 5, 4),
  B6 = c(5, 3, 3, 6, 5, 7, 7, 6, 6, 6),
  W6 = c(6, 5, 6, 7, 6, 5, 5, 7, 7, 5),
  B7 = c(2, 1, 3, 1, 3, 2, 3, 3, 2, 2),
  W7 = c(7, 7, 7, 7, 7, 1, 7, 1, 7, 7))
dat2 <- bws.dataset(respondent.dataset = res2,
                    response.type = 2,  # item number format
                    choice.sets = bibd,
                    design.type = 2,    # BIBD
                    item.names = items)
identical(dat, dat2) # check whether dat is identical to dat2


###########################################################
### 4.3.5 Analyzing responses using the counting approach
###########################################################
scores <- bws.count(data = dat)
scores

scores$aggregate$BW


###########################################################
### 4.3.6 Analyzing responses using the modeling approach
###########################################################
fr <- RES ~ A + B + C + E + F + G + strata(STR)
clogit(formula = fr, data = dat)


###########################################################
### 4.4.1 BWS based on a two-level OMED
###########################################################
sets.mfa <- cbind( # 12 rows x 9 columns
  c(1, 2, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2), # 1st column
  c(2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1), # 2nd column
  c(1, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 1), # 3rd column
  c(1, 2, 2, 2, 1, 2, 2, 1, 1, 1, 2, 1), # 4th column
  c(2, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2), # 5th column
  c(1, 1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1), # 6th column
  c(2, 1, 1, 2, 2, 2, 2, 1, 2, 1, 1, 1), # 7th column
  c(2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2), # 8th column
  c(2, 2, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1)) # 9th column
sets.mfa

items.mfa <- c(
  "Landscape",         # 1st column in sets.mfa
  "Biodiversity",      # 2nd column in sets.mfa
  "Water use",         # 3rd column in sets.mfa
  "Land conservation", # 4th column in sets.mfa
  "Flood control",     # 5th column in sets.mfa
  "Rural viability",   # 6th column in sets.mfa
  "Food security",     # 7th column in sets.mfa
  "Animal welfare",    # 8th column in sets.mfa
  "Cultural heritage") # 9th column in sets.mfa

bws.questionnaire(choice.sets = sets.mfa,
                  design.type = 1,
                  item.names = items.mfa)

data(mfa)
head(mfa) # display first 6 rows

data.mfa <- bws.dataset(
              respondent.dataset = mfa,
              response.type = 1, # row number format
              choice.sets = sets.mfa,
              design.type = 1)   # OMED
head(data.mfa)

count.mfa <- bws.count(data = data.mfa)
count.mfa

mar <- par("mar") # not shown in SPMUR
par(mar = c(5, 9, 2, 2)) # change margins of plot
o <- order(count.mfa$aggregate$stdBW) # permutation of stdBW
barplot(height = count.mfa$aggregate$stdBW[o],
        xlim = c(-1, 1),          # limits of x-axis
        names.arg = items.mfa[o], # names of each bar
        horiz = TRUE,   # bars are drawn horizontally
        las = 2)        # labels are perpendicular to axis
par(mar = mar) # not shown in SPMUR

fr.mfa <- RES ~ ITEM1 + ITEM2 + # exclude ITEM6 to
                ITEM3 + ITEM4 + #  normalize its
                ITEM5 + ITEM7 + #  coefficient to 0
                ITEM8 + ITEM9 +  strata(STR)
clg.mfa <- clogit(formula = fr.mfa, data = data.mfa)
clg.mfa

clg.mfa$coefficients
clgcoef.mfa <- c(clg.mfa$coef[1:5],
                 ITEM6 = 0,
                 clg.mfa$coef[6:8])
clgcoef.mfa
comp.mfa <- cbind(clogit = clgcoef.mfa, 
                  stdBW = count.mfa$aggregate$stdBW)
rownames(comp.mfa) <- items.mfa
comp.mfa
plot(comp.mfa, xlim = c(-1.5, 1.6))
text(comp.mfa, pos = 4,  # show names of roles to
     labels = items.mfa) #  the right of the points
cor(comp.mfa)

exp.clgcoef.mfa <- exp(clgcoef.mfa) # numerator
den <- sum(exp.clgcoef.mfa)         # denominator
sp <- matrix(exp.clgcoef.mfa / den, # shares of preference
             nrow = 9, ncol = 1,    # 9 rows x 1 column
             dimnames = list(items.mfa,    # set row and 
                             c("shares"))) #  column names
addmargins(sp, 1) # show shares of preferences with their sum


###########################################################
### 4.4.2 BWS based on a BIBD
###########################################################
sets.fruit <- cbind( # 7 rows x 4 columns
  c(1, 2, 2, 1, 1, 3, 1), # 1st column
  c(4, 3, 4, 2, 3, 5, 2), # 2nd column
  c(6, 4, 5, 5, 4, 6, 3), # 3rd column
  c(7, 6, 7, 6, 5, 7, 7)) # 4th column
sets.fruit

items.fruit <- c(
  "Apple",  # corresponds to value "1" in sets.fruit
  "Orange", # corresponds to value "2" in sets.fruit
  "Grapes", # corresponds to value "3" in sets.fruit
  "Banana", # corresponds to value "4" in sets.fruit
  "Peach",  # corresponds to value "5" in sets.fruit
  "Melon",  # corresponds to value "6" in sets.fruit
  "Pear")   # corresponds to value "7" in sets.fruit

bws.questionnaire(choice.sets = sets.fruit,
                  design.type = 2, # BIBD
                  item.names = items.fruit)

data(fruit)
head(fruit)

data.fruit <- bws.dataset(
                respondent.dataset = fruit,
                response.type = 1,  # row number format
                choice.sets = sets.fruit,
                design.type = 2,    # BIBD
                item.names = items.fruit)
head(data.fruit)

fr.fruit <- RES ~ Apple + Orange +  # exclude Pear to
                  Grapes + Banana + #  normalize its
                  Peach + Melon +   #  coefficient to 0
                  strata(STR)
clg.fruit <- clogit(formula = fr.fruit, data = data.fruit)
clg.fruit

count.fruit <- bws.count(data = data.fruit)
count.fruit

BW <- data.frame(count.fruit$disaggregate$BW)
BW <- lapply(BW, factor,            # scores are transformed
             levels = c(-4:4))      #  into a factor
score.tables <- lapply(BW, table)   # table for each item
ylimit <- max(unlist(score.tables)) # max count among items
par(mfrow = c(4,2)) # create 4 x 2 figure regions
for(i in 1:7){ # for each item
  barplot(height = score.tables[[i]],
          main = items.fruit[i], # title
          xlab = c("BW score"),  # label for x-axis
          ylab = c("Count"),     # label for y-axis
          ylim = c(0, ylimit))   # limits for y-axis
}
par(mfrow = c(1,1)) # not shown in SPMUR

Mean.BW <- colMeans(count.fruit$disaggregate$BW)
Standard.Deviation.BW <- apply(count.fruit$disaggregate$BW,
                               2, sd)
plot(x = Mean.BW,
     y = Standard.Deviation.BW, 
     xlim = c(-3, 3),
     ylim = c(1, 2))
text(x = Mean.BW,
     y = Standard.Deviation.BW,
     pos = 4,
     labels = items.fruit)

set.seed(123)
kmeans.fruit <- kmeans(x = count.fruit$disaggregate$BW,
                       centers = 2)
kmeans.fruit$centers
barplot(height = kmeans.fruit$centers, 
        names.arg = items.fruit, 
        beside = TRUE, # draw juxtaposed bar
        legend = c("cluster.1", "cluster.2"), 
        horiz = TRUE,
        las = 2,
        xlab = "BW score", 
        xlim = c(-4, 4))


### END ###