Stated Preference Methods Using R

###link code: http://lab.agr.hokudai.ac.jp/spmur/examples_ch2.R http://lab.agr.hokudai.ac.jp/spmur/

install.packages("DCchoice",
  repos = c("http://www.bioconductor.org/packages/release/bioc",
            "https://cran.rstudio.com/"),
  dep = TRUE)

## Installing package into 'C:/Users/binht/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)

## Warning: unable to access index for repository http://www.bioconductor.org/packages/release/bioc/bin/windows/contrib/4.0:
##   cannot open URL 'http://www.bioconductor.org/packages/release/bioc/bin/windows/contrib/4.0/PACKAGES'

## package 'DCchoice' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\binht\AppData\Local\Temp\Rtmp2JCx9T\downloaded_packages

library(DCchoice)

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

##   bid1 bid2 R1 R2 income work age female married
## 1  300  600  1  1   5000    0  24      0       0
## 2  300  150  0  0  40000    1  30      1       1
## 3  400  200  0  0   5000    0  70      0       0

data(CarsonSB)
CarsonSB

##    T1   Y   N
## 1  10 178  86
## 2  30 138 129
## 3  60 129 126
## 4 120  88 169

data(CarsonDB)
CarsonDB

##    T1  TU TL  yy yn ny  nn
## 1  10  30  5 119 59  8  78
## 2  30  60 10  69 69 31  98
## 3  60 120 30  54 75 25 101
## 4 120 250 60  35 53 30 139

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

##   bid R1
## 1  10  1
## 2  10  1
## 3  10  1
## 4  10  1
## 5  10  1
## 6  10  1

table(sb.data)

##      R1
## bid     0   1
##   10   86 178
##   30  129 138
##   60  126 129
##   120 169  88

data(KR)
head(KR, 3)

##   bid1 R1
## 1  100  1
## 2  100  1
## 3  100  1

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

##   bid1 bidh bidl answers age    sex income
## 1    6   18    3      yy   1 female      2
## 2   48  120   24      yn   2   male      1
## 3   48  120   24      yn   2 female      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)

##   bid1 bidh bidl answers age    sex income bid2 R1 R2 Lbid1 Lbid2
## 1    6   18    3      yy   1 female      2   18  1  1 1.792 2.890
## 2   48  120   24      yn   2   male      1  120  1  0 3.871 4.787
## 3   48  120   24      yn   2 female      3  120  1  0 3.871 4.787

0.0.1 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)

## 
## Call:
## sbchoice(formula = R1 ~ sex + age + income | Lbid1, data = NP)
## 
## Formula:
## R1 ~ sex + age + income | Lbid1
## 
## Coefficients: 
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.35184    0.64755   3.632 0.000281 ***
## sexfemale   -0.61025    0.25045  -2.437 0.014827 *  
## age         -0.37090    0.08546  -4.340 1.42e-05 ***
## income       0.26038    0.10566   2.464 0.013730 *  
## log(bid)    -0.46244    0.16212  -2.852 0.004338 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Distribution: log-logistic 
## Number of Obs.: 312  
## log-likelihood: -190 
## pseudo-R^2: 0.1142 , adjusted pseudo-R^2: 0.0909 
## LR statistic: 49.1 on 4 DF, p-value: 0.000 
## AIC: 390.578 , BIC: 409.293 
## 
## Iterations: 4  
## Convergence: TRUE 
## 
## WTP estimates:
##  Mean : Inf (because of |beta_bid| < 1)
##  Mean : 26.3 (truncated at the maximum bid)
##  Mean : 46.9 (truncated at the maximum bid with adjustment)
##  Median : 28.1

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

## 
## Call:
## sbchoice(formula = R1 ~ sex + age + income | Lbid1, data = NP, dist = "log-normal")
## 
## Formula:
## R1 ~ sex + age + income | Lbid1
## 
## Coefficients: 
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.43387    0.38527   3.722 0.000198 ***
## sexfemale   -0.36925    0.15125  -2.441 0.014633 *  
## age         -0.22497    0.05097  -4.414 1.02e-05 ***
## income       0.15219    0.06213   2.450 0.014302 *  
## log(bid)    -0.28076    0.09794  -2.867 0.004149 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Distribution: log-normal 
## Number of Obs.: 312  
## log-likelihood: -190 
## pseudo-R^2: 0.1134 , adjusted pseudo-R^2: 0.0901 
## LR statistic: 48.7 on 4 DF, p-value: 0.000 
## AIC: 390.932 , BIC: 409.647 
## 
## Iterations: 4  
## Convergence: TRUE 
## 
## WTP estimates:
##  Mean : 15573 
##  Mean : 26.2 (truncated at the maximum bid)
##  Mean : 46.5 (truncated at the maximum bid with adjustment)
##  Median : 27.4

names(sb.normal)

##  [1] "coefficients" "call"         "formula"      "glm.out"      "glm.null"    
##  [6] "distribution" "covariates"   "bid"          "nobs"         "yn"          
## [11] "data.name"    "terms"        "contrasts"    "data"         "xlevels"

sb.normal$coefficients

## (Intercept)   sexfemale         age      income    log(bid) 
##       1.434      -0.369      -0.225       0.152      -0.281

names(summary(sb.normal))

##  [1] "coefficients"      "call"              "formula"          
##  [4] "glm.out"           "glm.null"          "distribution"     
##  [7] "covariates"        "bid"               "nobs"             
## [10] "yn"                "data.name"         "terms"            
## [13] "contrasts"         "data"              "xlevels"          
## [16] "glm.summary"       "glm.null.summary"  "loglik"           
## [19] "loglik.null"       "medianWTP"         "meanWTP"          
## [22] "trunc.meanWTP"     "adj.trunc.meanWTP" "psdR2"            
## [25] "adjpsdR2"          "LR.test"           "AIC"

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

##    AIC    BIC 
## 390.93 409.65

summary(sb.normal)$glm.summary$coef

##             Estimate Std. Error z value   Pr(>|z|)
## (Intercept)  1.43387   0.385267  3.7217 1.9785e-04
## sexfemale   -0.36925   0.151250 -2.4413 1.4633e-02
## age         -0.22497   0.050973 -4.4135 1.0171e-05
## income       0.15219   0.062128  2.4496 1.4302e-02
## log(bid)    -0.28076   0.097942 -2.8666 4.1487e-03

0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

0.0.2 2.4.2 Estimating WTP with DBDC data

0.0.2.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

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)

## 
## Call:
## dbchoice(formula = R1 + R2 ~ female + age + married + log(income) | log(bid1) + 
##     log(bid2), data = AP)
## 
## Formula:
## R1 + R2 ~ female + age + married + log(income) | log(bid1) + log(bid2)
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)  3.866525   1.133178   3.412  0.000645 ***
## female       0.346016   0.149000   2.322  0.020219 *  
## age         -0.024648   0.004743  -5.197 < 2.2e-16 ***
## married     -0.356667   0.165999  -2.149  0.031666 *  
## log.income.  0.298830   0.108978   2.742  0.006105 ** 
## log(bid)    -1.271842   0.064474 -19.727 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Distribution: log-logistic 
## Number of Obs.: 721 
## Log-likelihood: -883.7191 
## 
## LR statistic: 53.173 on 4 DF, p-value: 0.000 
## AIC: 1779.4381 , BIC: 1806.9219 
## 
## Iterations: 44 13 
## Convergence: TRUE 
## 
## WTP estimates:
##  Mean : 371.53 
##  Mean : 181.57 (truncated at the maximum bid)
##  Mean : 193.42 (truncated at the maximum bid with adjustment)
##  Median: 93.577

summary(db.normal)

## 
## Call:
## dbchoice(formula = R1 + R2 ~ female + age + married + log(income) | log(bid1) + 
##     log(bid2), data = AP, dist = "log-normal", par = db.logit$coefficients)
## 
## Formula:
## R1 + R2 ~ female + age + married + log(income) | log(bid1) + log(bid2)
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)  2.304996   0.661228   3.486  0.000490 ***
## female       0.190012   0.087879   2.162  0.030602 *  
## age         -0.013967   0.002714  -5.146 < 2.2e-16 ***
## married     -0.225372   0.098014  -2.299  0.021483 *  
## log.income.  0.175737   0.063670   2.760  0.005778 ** 
## log(bid)    -0.759193   0.036193 -20.976 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Distribution: log-normal 
## Number of Obs.: 721 
## Log-likelihood: -879.4200 
## 
## LR statistic: 49.693 on 4 DF, p-value: 0.000 
## AIC: 1770.8400 , BIC: 1798.3238 
## 
## Iterations: 71 11 
## Convergence: TRUE 
## 
## WTP estimates:
##  Mean : 219.28 
##  Mean : 177.8 (truncated at the maximum bid)
##  Mean : 187.23 (truncated at the maximum bid with adjustment)
##  Median: 92.1

names(db.logit)

##  [1] "f.stage"      "dbchoice"     "coefficients" "call"         "formula"     
##  [6] "Hessian"      "distribution" "loglik"       "convergence"  "niter"       
## [11] "nobs"         "covariates"   "bid"          "yn"           "data.name"   
## [16] "terms"        "contrasts"    "data"         "xlevels"

coefficients(db.logit)

## (Intercept)      female         age     married log.income.    log(bid) 
##    3.866525    0.346016   -0.024648   -0.356667    0.298830   -1.271842

db.logit$loglik

## [1] -883.72

sum_db.logit <- summary(db.logit)
names(sum_db.logit)

##  [1] "f.stage"           "dbchoice"          "coefficients"     
##  [4] "call"              "formula"           "Hessian"          
##  [7] "distribution"      "loglik"            "convergence"      
## [10] "niter"             "nobs"              "covariates"       
## [13] "bid"               "yn"                "data.name"        
## [16] "terms"             "contrasts"         "data"             
## [19] "xlevels"           "medianWTP"         "meanWTP"          
## [22] "trunc.meanWTP"     "adj.trunc.meanWTP" "LR.test"          
## [25] "coef"              "AIC"

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

##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   3.8665    1.13318    3.41 0.000645
## female        0.3460    0.14900    2.32 0.020219
## age          -0.0246    0.00474   -5.20 0.000000
## married      -0.3567    0.16600   -2.15 0.031666
## log.income.   0.2988    0.10898    2.74 0.006105
## log(bid)     -1.2718    0.06447  -19.73 0.000000

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)

0.0.2.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.2

0.0.3 2.4.3 Constructing confidence intervals

0.0.3.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

options(digits = 3) # not shown in SPMUR

sblogit.krCI <- krCI(sb.logit)
sblogit.krCI

## the Krinsky and Robb simulated confidence intervals
##                         Estimate       LB       UB
## Mean                         Inf -999.000 -999.000
## truncated Mean            26.327   23.588   29.037
## adjusted truncated Mean   46.892   37.166   61.587
## Median                    28.138   16.175  113.475

sblogit.boCI <- bootCI(sb.logit)
sblogit.boCI

## the bootstrap confidence intervals
##                         Estimate       LB       UB
## Mean                         Inf -999.000 -999.000
## truncated Mean            26.327   23.346   29.432
## adjusted truncated Mean   46.892   36.124   63.895
## Median                    28.138   15.889  101.440

names(sblogit.krCI)

## [1] "mWTP"        "tr.mWTP"     "adj.tr.mWTP" "medWTP"      "out"

names(sblogit.boCI)

## [1] "mWTP"        "tr.mWTP"     "adj.tr.mWTP" "medWTP"      "out"

mean(sblogit.boCI$tr.mWTP)

## [1] 26.4

sd(sblogit.boCI$tr.mWTP)

## [1] 1.59

hist(sblogit.krCI$tr.mWTP, main="", 
     xlab="truncated mean WTP")

0.0.3.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.2

0.0.4 2.5.1 Kristrom’s nonparametric estimation of SBDC data

0.0.4.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

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

## Survival probability: 
##    Upper   Prob.
## 1      0  1.0000
## 2    100  0.4011
## 3    200  0.4011
## 4    300  0.3256
## 5    400  0.2294
## 6    Inf  0.0000
## 
## WTP estimates:
##  Mean: 135.71040  (Kaplan-Meier)
##  Mean: 201.60279  (Spearman-Karber)
##  Median: 83.48018

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

##      yes total yes/total
## 100   51    60     0.850
## 400   29    52     0.558
## 700   33    60     0.550
## 1000  31    57     0.544
## 1500  25    64     0.391
## 2000  16    56     0.286
## 2500  21    53     0.396
## 3000  16    53     0.302
## 5000  21    62     0.339
## 7000   5    45     0.111

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

## Survival probability: 
##     Upper   Prob.
## 1       0  1.0000
## 2     100  0.8500
## 3     400  0.5577
## 4     700  0.5500
## 5    1000  0.5439
## 6    1500  0.3906
## 7    2000  0.3394
## 8    2500  0.3394
## 9    3000  0.3217
## 10   5000  0.3217
## 11   7000  0.1111
## 12    Inf  0.0000
## 
## WTP estimates:
##  Mean: 2141.79768  (Kaplan-Meier)
##  Mean: 2519.99054  (Spearman-Karber)
##  Median: 1143.11270

plot(kr.example)

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

## 
## Call:
## sbchoice(formula = R1 ~ 1 | log(bid1), data = KR)
## 
## Formula:
## R1 ~ 1 | log(bid1)
## 
## Coefficients: 
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.63716    0.61741   7.511 5.88e-14 ***
## log(bid)    -0.68003    0.08493  -8.007 1.17e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Distribution: log-logistic 
## Number of Obs.: 562  
## log-likelihood: -346.936 
## pseudo-R^2: 0.1004 , adjusted pseudo-R^2: 0.0952 
## LR statistic: 77.456 on 1 DF, p-value: 0.000 
## AIC: 697.87247 , BIC: 706.53547 
## 
## Iterations: 4  
## Convergence: TRUE 
## 
## WTP estimates:
##  Mean : Inf (because of |beta_bid| < 1)
##  Mean : 2367.43 (truncated at the maximum bid)
##  Mean : 2960.86 (truncated at the maximum bid with adjustment)
##  Median : 915.076

0.0.4.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.2

0.0.5 2.5.2 Kaplan-Meier-Turnbull estimation of SBDC data

0.0.5.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

tb.NP <- turnbull.sb(R1 ~ bid1, data = NP)
summary(tb.NP)

## Survival probability: 
##    Upper   Prob.
## 1      0  1.0000
## 2      6  0.6579
## 3     12  0.5584
## 4     24  0.5122
## 5     48  0.4675
## 6    Inf  0.0000
## 
## WTP estimates:
##  Mean: 24.66514  (Kaplan-Meier)
##  Mean: 26.80324  (Spearman-Karber)
##  Median in: [     24 ,     48 ]

plot(tb.NP)

carson.sb <- turnbull.sb(R1 ~ bid, data = sb.data)
summary(carson.sb)

## Survival probability: 
##    Upper   Prob.
## 1      0  1.0000
## 2     10  0.6742
## 3     30  0.5169
## 4     60  0.5059
## 5    120  0.3424
## 6    Inf  0.0000
## 
## WTP estimates:
##  Mean: 52.80072  (Kaplan-Meier)
##  Mean: 61.07206  (Spearman-Karber)
##  Median in: [     60 ,    120 ]

plot(carson.sb) # not shown in SPMUR

0.0.5.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.2

0.0.6 2.5.3 Kaplan-Meier-Turnbull estimation of DBDC data

0.0.6.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

names(NP)

##  [1] "bid1"    "bidh"    "bidl"    "answers" "age"     "sex"     "income" 
##  [8] "bid2"    "R1"      "R2"      "Lbid1"   "Lbid2"

trdb.NP <- turnbull.db(R1 + R2 ~ bid1 + bid2, data = NP)
summary(trdb.NP)

## Survival probability: 
##    Upper    Prob.
## 1      0  1.00000
## 2      3  0.69831
## 3      6  0.65896
## 4     12  0.59574
## 5     18  0.41667
## 6     24  0.41667
## 7     48  0.22549
## 8    120  0.01879
## 9    Inf  0.00000
## 
## WTP estimates:
##  Mean: 19.41095  (Kaplan-Meier)
##  Mean: 30.38469  (Spearman-Karber)
##  Median in: [     12 ,     18 ]

plot(trdb.NP) # not shown in SPMUR