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


###########################################################
### 2.3.2 Installing the DCchoice package
###########################################################
# install.packages("DCchoice", 
#  repos = c("http://R-Forge.R-project.org", 
#            "http://www.bioconductor.org/packages/release/bioc",
#            "http://cran.rstudio.com/"), 
#  dep = TRUE)


###########################################################
### 2.3.3 Loading the package
###########################################################
library(DCchoice)    # loading the package


###########################################################
### 2.3.5 Example data sets
###########################################################
data(AP)     # loading the Albemarle-Pamlico sounds CV data
head(AP, 3)  # showing the first three rows

data(CarsonSB)
CarsonSB

data(CarsonDB)
CarsonDB

n <- rowSums(CarsonSB[, -1])
sb.data <- data.frame(
  bid = c(rep(CarsonSB$T1[1], n[1]), 
          rep(CarsonSB$T1[2], n[2]), 
          rep(CarsonSB$T1[3], n[3]), 
          rep(CarsonSB$T1[4], n[4])), 
  R1 = c(rep(1, CarsonSB$Y[1]), rep(0, CarsonSB$N[1]),
         rep(1, CarsonSB$Y[2]), rep(0, CarsonSB$N[2]),
         rep(1, CarsonSB$Y[3]), rep(0, CarsonSB$N[3]),
         rep(1, CarsonSB$Y[4]), rep(0, CarsonSB$N[4]))
)

head(sb.data)    # showing the first six rows
table(sb.data)

data(KR)
head(KR, 3)

data(NaturalPark, package = "Ecdat")
head(NaturalPark, 3)

NP <- NaturalPark	# a new object for transformation
NP$bid2 <- ifelse(NP$answers == "yy" | NP$answers == "yn", 
    NP$bidh, NP$bidl)
NP$R1 <- ifelse(NP$answers == "yy" | NP$answers == "yn",
                1, 0)
NP$R2 <- ifelse(NP$answers == "yy" | NP$answers == "ny",
                1, 0)
NP$Lbid1 <- log(NP$bid1)
NP$Lbid2 <- log(NP$bid2)

options(digits = 4) # not shown in SPMUR
head(NP, 3)


###########################################################
### 2.4.1 Estimating WTP with SBDC data
###########################################################
AP$Lbid1 <- log(AP$bid1)  # not shown in SPMUR
AP$Lbid2 <- log(AP$bid2)  # not shown in SPMUR
options(digits = max(3, getOption("digits") - 1)) # not shown in SPMUR

sb.logit <- sbchoice(R1 ~ sex + age + income | Lbid1, 
                     data = NP)
summary(sb.logit)

sb.normal <- sbchoice(R1 ~ sex + age + income | Lbid1, 
                      data = NP, 
                      dist = "log-normal")
summary(sb.normal)
names(sb.normal)
sb.normal$coefficients
names(summary(sb.normal))

options(digits = 5) # not shown in SPMUR
summary(sb.normal)$AIC
summary(sb.normal)$glm.summary$coef


###########################################################
### 2.4.2 Estimating WTP with DBDC data
###########################################################
data(AP)
db.logit <- dbchoice(R1 + R2 ~ female + age + married + 
                     log(income) | log(bid1) + log(bid2), 
                     data = AP)
db.normal <- dbchoice(R1 + R2 ~ female + age + married + 
                      log(income) | log(bid1) + log(bid2), 
                      data = AP, 
                      dist = "log-normal", 
                      par = db.logit$coefficients)
summary(db.logit)
summary(db.normal)
names(db.logit)
coefficients(db.logit)
db.logit$loglik
sum_db.logit <- summary(db.logit)
names(sum_db.logit)

options(digits = 3) # not shown in SPMUR
sum_db.logit$coef

options(digits = max(3, getOption("digits") - 1)) # not shown in SPMUR
write.table(sum_db.logit$coef, file = "dbout.txt", 
    quote = FALSE)
db.logit.null <- dbchoice(R1 + R2 ~ 1 | log(bid1) + 
                          log(bid2), data = AP)


###########################################################
### 2.4.3 Constructing confidence intervals
###########################################################
options(digits = 3) # not shown in SPMUR

sblogit.krCI <- krCI(sb.logit)
sblogit.krCI
sblogit.boCI <- bootCI(sb.logit)
sblogit.boCI
names(sblogit.krCI)
names(sblogit.boCI)
mean(sblogit.boCI$tr.mWTP)
sd(sblogit.boCI$tr.mWTP)
hist(sblogit.krCI$tr.mWTP, main="", 
     xlab="truncated mean WTP")


###########################################################
### 2.5.1 Kristrom's nonparametric estimation of SBDC data
###########################################################
options(digits = 6) # not shown in SPMUR
kr.AP <- kristrom(R1 ~ bid1, data = AP)
summary(kr.AP)
plot(kr.AP)

data(KR)
tab <- table(KR)
prop.tab <- prop.table(tab, margin = 1)
tab <- addmargins(tab, margin = 2)
tab <- cbind(tab[, -1], round(prop.tab[, 2], 3))
colnames(tab) <- c("yes", "total", "yes/total")
tab

options(digits = 6) # not shown in SPMUR
kr.example <- kristrom(R1 ~ bid1, data = KR)
summary(kr.example)

plot(kr.example)

kr.sb <- sbchoice(R1 ~ 1 | log(bid1), data = KR)
summary(kr.sb)


###########################################################
### 2.5.2 Kaplan-Meier-Turnbull estimation of SBDC data
###########################################################
tb.NP <- turnbull.sb(R1 ~ bid1, data = NP)
summary(tb.NP)
plot(tb.NP)

carson.sb <- turnbull.sb(R1 ~ bid, data = sb.data)
summary(carson.sb)
plot(carson.sb) # not shown in SPMUR


###########################################################
### 2.5.3 Kaplan-Meier-Turnbull estimation of DBDC data
###########################################################
names(NP)
trdb.NP <- turnbull.db(R1 + R2 ~ bid1 + bid2, data = NP)
summary(trdb.NP)
plot(trdb.NP) # not shown in SPMUR


###########################################################
### 2.A Appendix
###########################################################
nobs <- sum(CarsonDB[, 4:7])
db.data <- data.frame(bid1 = numeric(nobs),
    bid2 = numeric(nobs),
    R1 = numeric(nobs), R2 = numeric(nobs),
    bid2U = numeric(nobs), bid2L = numeric(nobs))
n <- rowSums(CarsonDB[, 4:7])
id2 <- cumsum(n)
id1 <- c(1, id2[1:3] + 1)

db.data$bid1[id1[1]:id2[1]] <- CarsonDB[1, 1]
db.data$bid1[id1[2]:id2[2]] <- CarsonDB[2, 1]
db.data$bid1[id1[3]:id2[3]] <- CarsonDB[3, 1]
db.data$bid1[id1[4]:id2[4]] <- CarsonDB[4, 1]
db.data$bid2U[id1[1]:id2[1]] <- CarsonDB[1, 2]
db.data$bid2U[id1[2]:id2[2]] <- CarsonDB[2, 2]
db.data$bid2U[id1[3]:id2[3]] <- CarsonDB[3, 2]
db.data$bid2U[id1[4]:id2[4]] <- CarsonDB[4, 2]
db.data$bid2L[id1[1]:id2[1]] <- CarsonDB[1, 3]
db.data$bid2L[id1[2]:id2[2]] <- CarsonDB[2, 3]
db.data$bid2L[id1[3]:id2[3]] <- CarsonDB[3, 3]
db.data$bid2L[id1[4]:id2[4]] <- CarsonDB[4, 3]

y1 <- CarsonDB$yy[CarsonDB$T1 == 10] +
      CarsonDB$yn[CarsonDB$T1 == 10]
y2 <- CarsonDB$yy[CarsonDB$T1 == 30] +
      CarsonDB$yn[CarsonDB$T1 == 30]
y3 <- CarsonDB$yy[CarsonDB$T1 == 60] +
      CarsonDB$yn[CarsonDB$T1 == 60]
y4 <- CarsonDB$yy[CarsonDB$T1 == 120] +
      CarsonDB$yn[CarsonDB$T1 == 120]
z1 <- CarsonDB$yy[CarsonDB$T1 == 10]
z2 <- CarsonDB$yy[CarsonDB$T1 == 30]
z3 <- CarsonDB$yy[CarsonDB$T1 == 60]
z4 <- CarsonDB$yy[CarsonDB$T1 == 120]
w1 <- CarsonDB$ny[CarsonDB$T1 == 10]
w2 <- CarsonDB$ny[CarsonDB$T1 == 30]
w3 <- CarsonDB$ny[CarsonDB$T1 == 60]
w4 <- CarsonDB$ny[CarsonDB$T1 == 120]

db.data$R1[db.data$bid1 == 10][1:y1] <- 1
db.data$R1[db.data$bid1 == 30][1:y2] <- 1
db.data$R1[db.data$bid1 == 60][1:y3] <- 1
db.data$R1[db.data$bid1 == 120][1:y4] <- 1
db.data$R2[db.data$R1 == 1 & db.data$bid1 == 10][1:z1] <- 1
db.data$R2[db.data$R1 == 1 & db.data$bid1 == 30][1:z2] <- 1
db.data$R2[db.data$R1 == 1 & db.data$bid1 == 60][1:z3] <- 1
db.data$R2[db.data$R1 == 1 & db.data$bid1 == 120][1:z4] <- 1
db.data$R2[db.data$R1 == 0 & db.data$bid1 == 10][1:w1] <- 1
db.data$R2[db.data$R1 == 0 & db.data$bid1 == 30][1:w2] <- 1
db.data$R2[db.data$R1 == 0 & db.data$bid1 == 60][1:w3] <- 1
db.data$R2[db.data$R1 == 0 & db.data$bid1 == 120][1:w4] <- 1
db.data$bid2[db.data$R1 == 1] <- db.data$bid2U[db.data$R1 == 1]
db.data$bid2[db.data$R1 == 0] <- db.data$bid2L[db.data$R1 == 0]


###########################################################
### 2.5.3 Kaplan-Meier-Turnbull estimation of DBDC data
###########################################################
head(db.data, 3)
carson.db <- turnbull.db(R1 + R2 ~ bid1 + bid2, 
                         data = db.data)
summary(carson.db)
plot(carson.db)


### END ###