Hi, and welcome!
A reproducible example, called a reprex is needed to provide a more useful answer. Without knowing the data that the function operates on, possible problems aren't detectable.
However, that may not be possible in this case. The DCchoice package is at version 0.0.15, relies on the version package that has an archived dependency on Icens, which can be installed via BiocManager::install("Icens") after installing BiocManager. If that weren't enough, the help(dbchoice) example requires installing the Ecdat package.
After those obstacles, considerable tinkering is required to produce a reprex from the example
library(DCchoice)
#> Loading required package: MASS
#> Loading required package: interval
#> Loading required package: survival
#> Loading required package: perm
#> Loading required package: Icens
#> Loading required package: MLEcens
#> Loading required package: Formula
require("Ecdat")
#> Loading required package: Ecdat
#> Loading required package: Ecfun
#>
#> Attaching package: 'Ecfun'
#> The following object is masked from 'package:base':
#>
#> sign
#>
#> Attaching package: 'Ecdat'
#> The following object is masked from 'package:MASS':
#>
#> SP500
#> The following object is masked from 'package:datasets':
#>
#> Orange
# may be necessary to install missing dependency of `intervals` dependency
# with BiocManager::install("Icens")
## Examples are based on a data set NaturalPark in the package
## Ecdat (Croissant 2011): DBDCCV style question for measuring
## willingness to pay for the preservation of the Alentejo Natural
## Park. The data set (dataframe) contains seven variables:
## bid1 (bid in the initial question), bidh (higher bid in the follow-up
## question), bidl (lower bid in the follow-up question), answers
## (response outcomes in a factor format with four levels of "nn",
## "ny", "yn", "yy"), respondents' characteristic variables such
## as age, sex and income (see NaturalPark for details).
data(NaturalPark, package = "Ecdat")
head(NaturalPark)
#> 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
#> 4 24 48 12 nn 5 female 1
#> 5 24 48 12 ny 6 female 2
#> 6 12 24 6 nn 4 male 2
## The variable answers are converted into a format that is suitable for the
## function dbchoice() as follows:
NaturalPark$R1 <- ifelse(substr(NaturalPark$answers, 1, 1) == "y", 1, 0)
NaturalPark$R2 <- ifelse(substr(NaturalPark$answers, 2, 2) == "y", 1, 0)
## We assume that the error distribution in the model is a
## log-logistic; therefore, the bid variables bid1 is converted
## into LBD1 as follows:
NaturalPark$LBD1 <- log(NaturalPark$bid1)
## Further, the variables bidh and bidl are integrated into one
## variable (bid2) and the variable is converted into LBD2 as follows:
NaturalPark$bid2 <- ifelse(NaturalPark$R1 == 1, NaturalPark$bidh, NaturalPark$bidl)
NaturalPark$LBD2 <- log(NaturalPark$bid2)
## The utility difference function is assumed to contain covariates (sex, age, and
## income) as well as two bid variables (LBD1 and LBD2) as follows:
fmdb <- R1 + R2 ~ sex + age + income | LBD1 + LBD2
## Not run:
## The formula may be alternatively defined as
## fmdb <- R1 + R2 ~ sex + age + income | log(bid1) + log(bid2)
## End(Not run)
## The function dbchoice() with the function fmdb and the dataframe
## NP is executed as follows:
NPdb <- dbchoice(fmdb, data = NaturalPark)
NPdb
#>
#> Distribution: log-logistic
#> (Intercept) sexfemale age income log(bid)
#> 3.490541 -0.267775 -0.351578 0.277374 -1.133728
NPdbs <- summary(NPdb)
NPdbs
#>
#> Call:
#> dbchoice(formula = fmdb, data = NaturalPark)
#>
#> Formula:
#> R1 + R2 ~ sex + age + income | LBD1 + LBD2
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.49054 0.45664 7.644 < 2.2e-16 ***
#> sexfemale -0.26778 0.21709 -1.234 0.217402
#> age -0.35158 0.07691 -4.571 5e-06 ***
#> income 0.27737 0.08747 3.171 0.001518 **
#> log(bid) -1.13373 0.08348 -13.581 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Distribution: log-logistic
#> Number of Obs.: 312
#> Log-likelihood: -398.892319
#>
#> LR statistic: 45.137 on 3 DF, p-value: 0.000
#> AIC: 807.784638 , BIC: 826.499654
#>
#> Iterations: 35 11
#> Convergence: TRUE
#>
#> WTP estimates:
#> Mean : 105.4603
#> Mean : 28.9656 (truncated at the maximum bid)
#> Mean : 31.45642 (truncated at the maximum bid with adjustment)
#> Median: 13.78239
## The confidence intervals for these WTPs are calculated using the
## function krCI() or bootCI() as follows:
## Not run:
krCI(NPdb)
#> the Krinsky and Robb simulated confidence intervals
#> Estimate LB UB
#> Mean 105.460 -484.020 1043.017
#> truncated Mean 28.966 25.313 33.583
#> adjusted truncated Mean 31.456 27.006 37.920
#> Median 13.782 11.358 16.617
# bootCI(NPdb) # stream of partial match warnings
## End(Not run)
## The WTP of a female with age = 5 and income = 3 is calculated
## using function krCI() or bootCI() as follows:
## Not run:
krCI(NPdb, individual = data.frame(sex = "female", age = 5, income = 3))
#> the Krinsky and Robb simulated confidence intervals
#> Estimate LB UB
#> Mean 58.0325 -172.7638 464.815
#> truncated Mean 19.0070 14.0477 25.282
#> adjusted truncated Mean 19.8373 14.4186 27.230
#> Median 7.5842 5.1206 11.148
# bootCI(NPdb, individual = data.frame(sex = "female", age = 5, income = 3))
## End(Not run)
## The variable age and income are deleted from the fitted model,
## and the updated model is fitted as follows:
update(NPdb, .~. - age - income |.)
#>
#> Distribution: log-logistic
#> (Intercept) sexfemale log(bid)
#> 2.921034 -0.403133 -1.025678
## The bid design used in this example is created as follows:
bid.design <- unique(NaturalPark[, c(1:3)])
bid.design <- log(bid.design)
colnames(bid.design) <- c("LBD1", "LBDH", "LBDL")
bid.design
#> LBD1 LBDH LBDL
#> 1 1.791759 2.890372 1.098612
#> 2 3.871201 4.787492 3.178054
#> 4 3.178054 3.871201 2.484907
#> 6 2.484907 3.178054 1.791759
## Respondents' utility and probability of choosing Yes-Yes, Yes-No,
## No-Yes, and No-No under the fitted model and original data are
## predicted as follows:
head(predict(NPdb, type = "utility", bid = bid.design))
#> f u l
#> [1,] 1.3945701 0.1490430 2.1804101
#> [2,] -1.3241270 -2.3629511 -0.5382869
#> [3,] -1.0371534 -2.0759774 -0.2513133
#> [4,] -1.8607951 -2.6466352 -1.0749550
#> [5,] -1.9349984 -2.7208384 -1.1491583
#> [6,] -0.1782278 -0.9640679 0.6076123
head(predict(NPdb, type = "probability", bid = bid.design))
#> yy nn yn ny
#> [1,] 0.53719194 0.1015235 0.2641289 0.09715566
#> [2,] 0.08604184 0.6314138 0.1240906 0.15845369
#> [3,] 0.11145371 0.5624997 0.1502459 0.17580065
#> [4,] 0.06619670 0.7455381 0.0684137 0.11985151
#> [5,] 0.06175487 0.7593571 0.0644435 0.11444449
#> [6,] 0.27606448 0.3526041 0.1794961 0.19183531
## Utility and probability of choosing Yes for a female with age = 5
## and income = 3 under bid = 10 are predicted as follows:
predict(NPdb, type = "utility",
newdata = data.frame(sex = "female", age = 5, income = 3, LBD1 = log(10)))
#> [1] -0.3135033
predict(NPdb, type = "probability",
newdata = data.frame(sex = "female", age = 5, income = 3, LBD1 = log(10)))
#> [1] 0.4222599
## Plot of probabilities of choosing yes is drawn as drawn as follows:
plot(NPdb)
## The range of bid can be limited (e.g., [log(10), log(20)]):
plot(NPdb, bid = c(log(10), log(20)))
Created on 2020-01-16 by the reprex package (v0.3.0)
So, we at least know that the function will work on the NaturalPark data.
That does not guarantee, however, that it will work on all data.
Looking under the hood with
dbchoice
we find that it's a wrapper on optim
dbchoice <- optim(ini, fn = dbLL, method = "BFGS",
hessian = TRUE, dvar = yvar, ivar = mmX, bid = bid, control = list(abstol = 10^(-30))))
This old S/O post explains the source of the problem, which is that when a zero value is passed to optim things fall apart.