Methods

To assess the effects of DCE design features and model assumptions on DCE results, we completed a series of computer simulations based on 2 actual DCEs. One of the co-investigators provided access to the data from the 2 DCEs, which we selected based on data availability and applicability to US health care. Study 1 examined preferences for organ allocation among adults (N = 2051) in the Australian general public.³⁷ That study used an efficient design, with respondents randomly assigned to 5 blocks of 30 choice tasks. Each task included 2 alternatives; the alternatives had 15 attributes, with 2 to 6 levels per attribute. (Appendix A: Table A1 lists the study 1 attributes and variable names.) Study 2 examined preferences for labor induction among women (N = 362) who were participating in a randomized trial of labor induction alternatives (n = 260) or who were pregnant and volunteered only for the DCE (n = 102).³⁸ That study used an efficient fractional factorial design, with participants randomly assigned to 2 blocks of 25 choice tasks. Each task included 3 alternatives (outpatient, basic inpatient, and enhanced inpatient care for labor induction). The alternatives, based on the actual treatment arms of the randomized trial, had 6 attributes, with 1 to 4 levels per attribute. In study 2, some attributes were alternative specific (eg, relevant only for a particular type of care), and there was a fixed status quo option (basic inpatient care). (In a DCE, a status quo option, if present, is a fixed set of attributes that respondents compare with the other alternatives presented. The status quo option often represents usual care or the current state of affairs. Appendix A: Table A2 lists the study 2 attributes and variable names.)

The purpose of conducting simulations based on these 2 studies was to use realistic data to explore the effects on parameter estimates and their standard errors under violations of certain assumptions or combinations thereof: an unbiased sample, independence of utility from the set of attributes presented (eg, omitting an attribute does not affect the utility of other attributes), no ignoring of attributes, and no unmeasured interactions (which could be tested by including an interaction in the data-generating process but not in the analysis model). We explored these effects with samples of various sizes. We had expectations for the effects of individual simulation parameters. For example, selection bias in which sample selection is positively associated with a given parameter should bias the estimate of that parameter upward; omission of attributes that cause variance inflation or deflation because of respondent behavior should lead to inflated or deflated SEs, respectively; ignoring attributes should decrease their apparent importance and increase the apparent importance of other attributes; and an unmeasured interaction should result in biased main-effect estimates for the variables involved in the interaction. However, because of the complexity of the simulation setup, we did not specify hypotheses for the combined effects of simulation parameters. Rather, we aimed to describe those effects.

We originally considered generating our simulated populations using 4 actual DCEs, which we selected based on data availability. The final research protocol provided for 3 actual DCEs and, as described in the section on stakeholder engagement, we reduced the number of studies to 2 to expand other aspects of the project. Compared with the studies examined in our systematic review, study 1 was atypical in that it had a fairly large sample, focused on developing a priority-setting framework (as did 9% of the studies reviewed), and had relatively large numbers of choice tasks and attributes. Study 2 had a more typical sample size, consumer experience focus, and number of attributes but still had a higher-than-average number of choice tasks. Of the studies that we dropped, both had small sample sizes (130 and 219). One had a common number of attributes (4) but many choice tasks (32); the other had many attributes (9) and choice tasks (32). The implications of our choice of studies are discussed in the Limitations section toward the end of the report.

Beginning with the actual data from each study, we verified that we could replicate the original results. We then simulated a population of 10 000 individuals based on the estimates from the original study. For each fixed parameter, the entire sample received the estimated mean parameter value. For each random parameter, we used the estimated mean and SD to simulate a normally distributed parameter across the population. The simulated population included an additional dichotomous covariate, x, which had 50% prevalence and was correlated with selected random parameters at 0, 0.20, or 0.50. We used the BinNor package in R to simultaneously induce correlations between x and those parameters.

In theory, high task complexity may cause respondents to ignore less important attributes, which would violate the assumptions of the usual DCE analysis models. To flag complex tasks, we used the parameter estimates from the original study to measure the entropy (complexity) of each task. Entropy measures the amount of information contained in a task. For a given choice task, entropy is calculated as $-{\displaystyle {\sum }_k{p}_k}\cdot \mathit{ln}(p_k)$, where p_k is the probability of choosing alternative k.³⁹ Tasks with different alternatives have smaller entropy values, whereas tasks with similar alternatives have larger entropy values. For example, in a 2-alternative task, entropy is maximized (at 0.69) when each alternative has probability 0.5 of being selected; if the probabilities are close to 0 and 1, as in a task with a dominant alternative, the entropy is near 0. Tasks with higher entropy (ie, tasks without a clearly dominant alternative) require more thinking. Therefore, we flagged tasks with entropy above the 75th percentile as complex. In our simulated populations, described further in the next paragraph, some people responded to these complex tasks by ignoring certain attributes.

For each study, we used North Carolina State University's high-performance computing cluster to simulate the usual DCE study process: a pilot study, followed by a redesign based on pilot parameter estimates, followed by a larger main DCE. We used the set of simulation parameters shown in Table 1. Starting with either study 1 or study 2 as a data source, we selected pilot samples of 4 sizes from the population. Pilot sample selection was either completely random or minimally or moderately biased (ie, correlated with the unmeasured variable x at r = 0.2 or r = 0.5, where x was correlated with selected preference parameters). Then, using the design of the real DCE, we simulated respondents' choices. We either (1) omitted selected attributes from the data-generating process and analysis model, in which case variance was either increased or decreased by 10% (because simpler tasks could either improve measurement precision by decreasing fatigue or worsen precision by withholding information that participants need to make choices); or (2) included those attributes in the data-generating process and analysis model, with or without an unmeasured interaction (because failing to account for an interaction could bias DCE results). When including the selected attributes in the data-generating process (condition 2 in the previous sentence), we set their parameter values such that the utility ranges associated with those attributes would be 20% of the maximum utility range of any attribute in the original study (ie, in a sense, the selected attributes would be 20% as important as the most important attribute). We randomly selected 50% of the study population to ignore the least important attributes (as measured using the parameter estimates from the original study) in complex tasks (those with entropy above the 75th percentile).

Table 1

Simulation Parameters.

In the simulated population of 100 000 people for study 1, a total of 49 984 people (50%) were randomly assigned to ignore selected attributes in high-entropy tasks, and 50 166 (50%) were randomly assigned a value of 1 for unmeasured variable x, which had a correlation of 0.35 with qol_post and a correlation of −0.35 with qol_pre. The population was generated using parameters with the following means (SDs for random parameters only): asc1 0.097, age5 0.673, age15 0.573, age25 0.388, age55 −0.283, age70 −1.179, prev −0.149, adher −0.063, sex 0.08, donor 0.199, depend 0.356, donor*depend 0.365, wait 0.043 (0.052), le_pre −0.089 (0.076), le_post 0.06 (0.05), qol_pre −0.059 (0.059), qol_post 0.114 (0.108), comb1 0.002, comb2 −0.068, comb3 −0.096, comb4 −0.23, smok1 −0.263, smok2 −0.787, drink1 −0.095, drink2 −0.359, and obes −0.272.

For study 2, in the simulated population of 100 000 people, 49 917 people (50%) were randomly assigned to ignore selected attributes in high-entropy tasks, and 50 297 (50%) were randomly assigned a value of 1 for unmeasured variable x, which had a correlation of 0.35 with midw. The population was generated using normally distributed parameters with the following means (SDs for random parameters only): asc1 1.635, env_home 0.324, env_pvt_shrd 0.886, env_pvt_pvt 1.8578, pain_1dose 0.178, pain_mild 0.606, pain_mild_nd 0.572, chk_ph_avl −0.791, chk_mw_avl 1.361, midw 0.173 (0.057), trip −0.265 (0.088), tt −0.014 (0.005), and trip*tt −0.005.

To increase the likelihood of stable estimation with 1000 Halton draws, we generated only 5 random parameters for each study. Because parameters had to vary to be correlated with unmeasured variable x (and therefore with sample selection in certain scenarios), we included as random parameters the variables that were correlated with x. After generating the simulated population, for each simulation scenario we estimated an MNL model (Appendix B: Table B1) with utility as a function of the attributes of alternatives. We conducted the pilot simulations using R version 3.5.1, with the gmnl package for MNL estimation. We used MNL estimation for the pilot simulations because the pilot sample sizes were insufficient to support estimation of RPL models. As an example, the pilot analysis model for study 1 in scenarios with no ignored attributes or omitted terms was:

where U_nk = utility for person n and alternative k; asc1 is an alternative-specific constant for alternative 1; n = 1, 2, …, N; and ε_nk is the independently and identically distributed random error, following the Gumbel distribution.

For the redesign step, we used SAS version 9.4 with a set of SAS macros to generate a new, efficient design for each simulation scenario based on the parameter estimates from each pilot study.⁴⁰ The redesign software was programmed to generate a D-efficient design for MNL estimation, assuming that the parameter estimates from the pilot study were correct and keeping constant the numbers of tasks, blocks, and alternatives.

We based the main DCE simulation on the new DCE design. In the main DCE simulation, we varied 2 parameters: sample size and analysis model. We randomly selected a sample of each size from the simulated population. The analysis models included a full model (main effects plus interaction), a model omitting the interaction, and a model omitting both the interaction and 1 of the main effects. For the main DCE, as for the pilot, we used the gmnl package in R version 3.5.1. Although we had originally planned to use both MNL and RPL models in the main DCE phase for comparison, we implemented only MNL models, which by definition were misspecified because they ignored preference heterogeneity. We made this simplification based on consultation with stakeholders, as described previously. As in the pilot phase, the MNL models expressed utility as a function of attributes of alternatives.

As an example, the full model for study 1 (main effects plus interaction) is as follows:

where U_nk = utility for person n and alternative k; asc1 is an alternative-specific constant for alternative 1; n = 1, 2, …, N; and ε_nk is the independently and identically distributed random error, following the Gumbel distribution.

Using study 1 as an example, Figure 2 shows the flow of the simulations, from simulation of the study population to pilot design to pilot simulation to main DCE simulation. In the pilot design phase, the parameters involving omission of selected attributes (with variance inflation or deflation), inclusion of those attributes (with or without an interaction in the data-generating process), and ignoring selected attributes in complex tasks account for 2 × 2 × 2 = 8 sets of simulation scenarios. In the pilot simulation phase, the sample size and sample selection simulation parameters account for 4 × 3 = 12 sets. Finally, in the main DCE phase, the sample size and analysis model simulation parameters account for 3 × 3 = 9 sets of simulation parameters. Therefore, the total number of simulation scenarios for each original study was 8 × 12 × 9 = 864. For each simulation scenario, we used 5000 iterations in the pilot simulation phase. Based on the results of each pilot study, we obtained a design for the main DCE phase in which we used 1000 iterations for study 1 and 5000 iterations for study 2. We used convergence plots to select the number of iterations, enabling us to conserve computing resources while ensuring that we had used enough iterations for stable estimation.

Figure 2

Simulation Flow Chart for Study 1.

In the convergence plots, we plotted parameter estimates, ratios of parameter estimates, and standard errors by number of iterations, from 100 up to the maximum number of iterations (1000 or 5000) in increments of 100. We produced convergence plots for “best-case” and “worst-case” simulation scenarios. The best case included pilot N = 240, all attributes included in the pilot data-generating process and model, no unmeasured interaction in the pilot DCE, no selection bias, no attributes ignored, and the correct main DCE analysis model. The worst case included pilot N = 30, selected attributes omitted from the pilot data-generating process and model, variance inflated by 10%, selection bias (r = 0.5), selected attributes ignored in complex tasks, and an unmeasured interaction in the main DCE analysis model.

To assess the combined effects of the pilot and main DCE simulation parameters, we estimated bias, relative standard error (RSE, or the standard error divided by the parameter estimate), and D-error (the determinant of the asymptotic covariance matrix, a measure of overall error in parameter estimation). Because DCE parameter estimates are relative rather than absolute, we estimated bias in ratios of parameter estimates (relative to qol_post in study 1 and midw in study 2) rather than estimating bias directly. For pairs of attributes that appear only in linear terms in the analysis model, this ratio is known as a marginal rate of substitution (MRS).⁷ We expressed each ratio as a percentage of the denominator (percent bias). (For qol_post and midw, respectively, we used qol_pre and env_pvt_pvt as the denominator when estimating bias.) High RSE values indicate imprecise estimates.

We then plotted bias, RSE, and D-error for different combinations of pilot and main DCE simulation parameters, holding other simulation parameters constant at their default values. We assigned the default values with the goal of making the comparisons as informative as possible. The default pilot sample size was 240. With all attributes included in the pilot data-generating process, the default was to include no unmeasured interaction; with selected attributes excluded, the default variance adjustment was to decrease variance by 10%. By default, selection bias was absent, no one ignored attributes in complex tasks, the main DCE sample size was 1000, and the main DCE included the full-analysis model. When comparing across different main DCE analysis models, we included an unmeasured interaction in the data-generating process (and not in the pilot analysis model).

Results

We completed the simulations without errors, except that for study 2, a coding error caused some lines of output to be garbled at random because multiple processes were writing to the same output file. This error affected <1% of output (18 815 of 4 296 865 lines). Because of computing resource limitations, we were unable to complete as many iterations for study 1 as for study 2. However, the convergence plots (Appendix A: Figures A1-A4) showed little fluctuation in estimates with 1000 to 5000 iterations.

The following subsections report the findings as percent bias, RSE, and D-error, using selected figures as examples. Full sets of figures appear in Appendix A. Each figure in this section represents a combined effect of pilot DCE simulation parameters and main DCE simulation parameters. For example, Figure 3 shows percent bias in the estimate of the parameter age70 (relative to the qol_post parameter estimate) by pilot sample size and main DCE sample size.

Figure 3

Percent Bias in the age70 to qol_post Ratio by Selected Study 1 Simulation Parameters.

Percent Bias and Main DCE Sample Size

Little bias emerged in the default scenario. (In study 1, an exception was that estimates for 1 of the comorbidity variables, comb1, were strongly biased. However, the true population mean of that parameter was 0.002—near 0—which may explain the large differences in the estimated vs actual parameter ratios.) Larger main DCE sample sizes generally resulted in percent bias closer to 0. We observed no consistent combined effect of pilot sample size and pilot selection bias on the one hand and main DCE sample size on the other hand (eg, Figure 3). (For example, the leftmost column of Appendix A: Figure A5 shows generally low levels of bias that appear to fluctuate randomly across pilot study sample sizes; for a main DCE sample size of 400, the percentage of plots in which the absolute value of bias increased from 1 pilot study sample size to the next—30 to 100, 100 to 170, or 170 to 240—was 52%, 56%, and 40%, respectively.) In study 1, we saw no consistent combined effect of the other 3 pilot simulation parameters (ignoring attributes, unmeasured interaction, and omission of attributes plus variance inflation) and the main DCE sample size (Appendix A: Figure A5). However, in study 2, we observed that when 50% of simulated respondents ignored the professional availability attributes (chk_ph_avl and chk_mw_avl) in complex tasks, many estimates were slightly biased, usually upward (eg, Figure 4; see also Appendix A: Figure A6). Given that bias was measured in reference to a ratio of parameter estimates, the bias and its direction were consistent, with an apparent decrease in the magnitude of the ignored parameters' utility and proportional increases in the magnitude of other parameters' utility. A similar pattern emerged in study 2 when the pilot DCE involved an unmeasured interaction (which is known to cause bias in the coefficients of variables correlated with the interaction term) and also when selected attributes (trip, tt) were omitted from the pilot study and the variance was inflated by 10% (Appendix A: Figure A6). When the trip and tt attributes were omitted in the pilot study, positive bias appeared in the ratios of most parameters to midw, with the possible exception of trip and tt themselves, and negative bias appeared in the ratio of midw to env_pvt_pvt. Omission of the trip and tt attributes would have led to 0 (ie, uninformative) priors for those 2 parameters in the redesign stage. This should not have biased the priors for other parameters, but greater measurement error for trip and tt in the main DCE phase may have inflated the apparent importance of other parameters in the main DCE analysis model.

Figure 4

Percent Bias in the chk_ph_avl to midw Ratio by Selected Study 2 Simulation Parameters.

Percent Bias and Main DCE Analysis Model

In study 1, misspecified models resulted in bias in certain parameter estimates. Model 1, which omitted an interaction, yielded biased estimates for the 2 terms (donor and depend) that were involved in the interaction (Appendix A: Figure A7). Model 0, which omitted the interaction and also the main effect of depend, yielded even greater bias in the estimate for donor (Figure 5) as well as bias in the alternative-specific constant and some of the comorbidity variables (Appendix A: Figure A7). We observed no clear combined effect of pilot and main DCE simulation parameters (Appendix A: Figure A7).

Figure 5

Percent Bias in the donor to qol_post Ratio by Selected Study 1 Simulation Parameters.

In study 2, we observed similar effects of main DCE model on the parameter estimates for the terms involved in the interaction (Appendix A: Figure A8). We also observed that when the interaction term was omitted (model 1), other parameter estimates reflected some bias (Appendix A: Figure A8). When a main effect was also omitted (model 0), many parameter estimates reflected large bias. In 2 instances, pilot simulation parameters appeared to have a combined effect with main DCE simulation parameters. First, when the main DCE model omitted both the interaction term (trip*tt) and a main effect (tt) (model 0), the amount of selection bias in the pilot study appeared to affect the amount of bias, increasing the positive bias in tt but often reducing the amount of negative bias in other parameter estimates (relative to midw) (eg, Figure 6). Second, when the pilot data-generating process and model omitted certain attributes and the pilot variance was inflated by 10%, the main DCE parameter estimates tended to be more positively biased (relative to midw) (Appendix A: Figure A8).

Figure 6

Percent Bias in the chk_ph_avl to qol_post Ratio by Selected Study 2 Simulation Parameters.

RSE and Main DCE Sample Size

In both study 1 and study 2, main DCE sample size was negatively associated with RSE, and all plots showed nearly parallel lines, indicating that pilot and main DCE simulation parameters had no combined effect on RSE (Appendix A: Figures A9 and A10).

RSE and Main DCE Analysis Model

In study 1, the main DCE analysis model appeared not to have a consistent effect on RSE (Appendix A: Figure A11). The full model (model 2) resulted in a higher RSE for both terms involved in the interaction (depend and donor) (eg, Figure 7). We observed no consistent combined effect of pilot and main DCE simulation parameters. In study 2, when the analysis model omitted the interaction (model 1), we observed a greater RSE for trip (but not tt) (Figure 8) and no effect on most other parameters (Appendix A: Figure A12). With both the interaction and the main effect of trip omitted (model 0), we observed a greater RSE for some parameters but a lower RSE for others (eg, trip, pain_mild_nd, midw). We saw no apparent combined effect of pilot and main DCE simulation parameters.

Figure 7

RSE for the donor Parameter Estimate by Selected Study 1 Simulation Parameters.

Figure 8

RSE for the trip Parameter Estimate by Selected Study 2 Simulation Parameters.

D-Error

In both study 1 (Appendix A: Figures A13 and A15) and study 2 (Appendix A: Figures A14 and A16), when the main DCE had a small sample size or the main DCE model omitted the interaction, we observed greater overall error in parameter estimation, especially when the pilot sample size was small, the interaction term was omitted from the pilot data-generating process and analysis model, or selected attributes were omitted from the pilot data-generating process and model and variance were reduced (vs increased) by 10% (eg, Figure 9). Especially in study 2, omitting both the interaction and a main effect from the main DCE analysis model resulted in a greater increase in overall error when the pilot sample size was small or the pilot study had strong selection bias (Figure 10).

Figure 9

D-Error by Selected Study 1 Simulation Parameters.

Figure 10

D-Error by Selected Study 2 Simulation Parameters.

Method

The purpose of this analysis was to demonstrate how Halton sequences correlate with each other because (1) correlations among Halton sequences could lead to correlations among random parameters, which would violate the assumptions of the RPL model; and (2) correlations among Halton sequences could cause the sequences to do a poor job of covering the 0 to 1 interval, which in turn would limit the simulated values of random parameters during the estimation process, possibly causing errors in parameter estimation. With R version 3.5.0,⁴¹ we used the randtoolbox package to generate 10 000 draws from each of 50 Halton sequences. Because we used 50 sequences and the 50th prime is 229, we skipped the first 229 elements of each sequence to avoid excessive correlations between sequences.²³ To show correlations between Halton sequences and bivariate coverage of the 0 to 1 interval, we created miniature scatter plots between pairs of adjacent Halton sequences with 500, 1000, 5000, and 10 000 draws. We used these plots to demonstrate how correlation and coverage varied with the number of Halton sequences and the number of draws from each sequence. For selected scatter plots that showed evidence of correlation and incomplete coverage, we also generated surface plots to demonstrate departures from bivariate normality.

Results

Figure 11 shows scatter plots of adjacent Halton sequences with different numbers of draws. In each panel, plot 1 represents Halton sequences 1 and 2, plot 2 represents sequences 2 and 3, and so forth. With 500 draws, increasing correlation and decreasing coverage of the 0 to 1 interval became apparent with 12 to 14 sequences (panel a, plots 11 and 13). With 1000 draws, correlation and coverage issues began with 18 sequences (panel c, plot 17). With 5000 and 10 000 draws, issues appeared once the number of sequences reached 34 (panel e, plot 33) and 42 (panel f, plot 41) respectively. We observed that the 4 plots referenced here are also associated with departures from bivariate normality (not shown). To summarize, Halton sequences based on larger primes showed greater correlation and more coverage issues; using more draws mitigated these issues.

Figure 11

Scatter Plots of the (k + 1)st Halton Sequence Against kth Halton Sequence for k = {1.49} and Number of Draws R = {500, 1000, 5000, 10 000}.

Method

Because the purpose of Halton sequences in RPL estimation is to simulate random parameters with specific distributions, the next step after generating the sequences is to convert the Halton draws to Halton-normal draws (evenly spaced draws from the normal distribution, generated using Halton sequences). We effected this conversion by treating the Halton sequences described in the previous subsection as lists of normal quantiles. We examined univariate and multivariate normality because the RPL model assumes specific distributions for random parameters, and violations of this assumption for normally distributed random parameters could lead to error in parameter estimation.

To examine univariate normality, we conducted Shapiro-Wilk tests for each of the 50 Halton sequences, using 50 to 5000 draws in increments of 50. We created heat maps showing, for each combination of number of draws and number of random parameters (Halton sequences), the percentage of tests in which the null hypothesis was retained, the minimum P value, and the median P value.

To examine multivariate normality, we conducted Henze-Zirkler tests for each of the 50 Halton sequences, using between 50 and 10 000 draws in increments of 50. We created a heat map showing, for each combination of number of draws and number of random parameters, the percentage of tests in which the null hypothesis was retained. We also created a line plot showing the minimum number of draws required for the Henze-Zirkler P value to remain above thresholds of P < .25 and P < .05.

Results: Univariate Normality

The heat map in Figure 12, panel (a) shows the percentage of Shapiro-Wilk tests for which the null hypothesis of univariate normality was retained, by number of draws and number of random parameters. Panels (b) and (c) show the minimum and median P values, respectively. With 5 random parameters, the P values for some Shapiro-Wilk tests began to decrease (panel b). With 500 draws and 10 random parameters, or 1000 and 12 respectively, at least 1 random parameter departed from normality (panels a and b). With >10 random parameters, the number of draws required to maintain univariate normality increased sharply (panels a and b); with 10 to 15 random parameters, using 1000 draws did not prevent a substantial proportion of random parameters from being non-normally distributed (panel a). With 500 draws and 17 random parameters, or 1000 and 22 respectively, half of the random parameters departed from normality (panel c).

Figure 12

Percentage of Univariate Normality, Minimum Shapiro-Wilk P Value, and Median Shapiro-Wilk P Value by Number of Halton Draws and Number of Random Parameters.

Results: Multivariate Normality

Figure 13 shows that as the number of random parameters increased to 7 to 12 and beyond, the P value for the Henze-Zirkler test of multivariate normality decreased, meaning that increasing the number of random parameters increased the likelihood of departures from multivariate normality. The number of draws required to keep the P value above .25 (Figure 14, dashed line) was 500 with 10 random parameters and increased sharply as additional random parameters were introduced. Similarly, the number of draws required to keep the P value above .05 (Figure 14, solid line) was 4000 with 11 random parameters and increased sharply with the number of random parameters. By overlaying Figure 1 onto Figure 14, we observed that, among the 40 recently published health-related DCEs that described RPL results and reported number of draws, the 40% (16/40) with ≥10 random parameters used too few draws to prevent a noticeable departure from the multivariate normal distribution. (We did not assess the distributional assumptions in those studies, however.)

Figure 13

Henze-Zirkler P Value by Number of Halton Draws and Number of Random Parameters.

Figure 14

Minimum Number of Halton Draws Needed for Henze-Zirkler P Value to Remain Above Selected Cutoffs, by Number of Halton Draws and Number of Random Parameters.

Method

The purpose of this analysis was to assess, for different numbers of Halton sequences (random parameters to be estimated), how many Halton draws were required to keep the correlations among draws below a given level. As mentioned previously, excessive correlations among random parameters would violate the assumptions of the RPL model and could cause errors in parameter estimation. Using the same data as for the multivariate normality tests (in the previous subsection), we created a heat map showing, for each combination of number of draws and number of random parameters, the maximum Spearman correlation between any pair of Halton-normal draws. We superimposed a line plot showing the minimum number of draws required for the maximum correlation to remain below thresholds of r < 0.2 and r < 0.1.

Results

Figure 15 shows that as the number of draws decreased or the number of random parameters increased, the maximum Spearman correlation between pairs of Halton-normal variables increased, meaning that using too few draws (or increasing the number of random parameters without increasing the number of draws) led to violations of the independence assumption. Keeping the maximum correlation below 0.2 required about 250 draws with 10 to 15 random parameters, 500 draws with 17 random parameters, and ≥1000 draws with ≥22 random parameters. Obeying the more conservative threshold of 0.1 required 500 draws with 13 random parameters and 1000 draws with 17 random parameters. With either threshold, the required number of draws increased sharply as the number of random parameters increased beyond 15 or 20.

Figure 15

Minimum Number of Halton Draws Needed for Correlations Between Halton-Normal Draws to Remain Below Selected Cutoffs, by Number of Halton Draws and Number of Random Parameters.

Overlaying Figure 1 onto Figure 15 showed that, among the 75% (30/40) of studies with ≤12 random parameters, most used enough draws to keep correlations among normally distributed random parameters below 0.1. Among the 10 studies with ≥13 random parameters, 7 could be expected to have correlations >0.2 if they assumed normally distributed random parameters.

Method

To examine the effect of correlated random parameters on bias and variance estimation, we used the R version 3.5.0 software⁴¹ to simulate a simple DCE with 3 alternatives, three 3-level attributes, a fractional-factorial orthogonal design with 6 choice tasks, and 500 participants. The DCE design did not include an opt-out alternative. We generated the data with no alternative-specific constant, no effect of the ordering of the alternatives on preferences, and a random parameter for each attribute (b1, b2, and b3). We set the true means at 1/3, 2/3, and 1.0 and the SDs at 0.5.

We created 10 simulation scenarios, with all pairs of random parameters (b1 and b2, b2 and b3, and b1 and b3) correlated at r = 0.0 to 0.9 by steps of 0.1. For each scenario, we used the R version 3.5.1 software⁴¹ to simulate 1000 RPL models, each using 1000 pseudorandom draws with the specified pairwise correlations. We used the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm in the R maxLik package⁴² for estimation and a different random number seed for each of the 1000 iterations. We estimated the 3 random parameters as linear effects, along with their SDs and the intercept. Therefore, the model (Appendix B: Table B1) estimated 7 parameters with (3 − 1) × 6 = 12 df. Across all iterations within each scenario, we measured for each parameter the mean bias (estimate minus true effect), mean squared error (MSE) (mean of squared bias), SE (SD of estimates), and 95% CI coverage (percentage of confidence intervals that contained the true value). We compared these statistics across correlation levels.

Results

For each level of simulated correlation among pairs of pseudorandom sequences, Figure 16 shows the mean bias, mean absolute bias, MSE, and standard error of each parameter estimate across simulation runs. The top and bottom rows of Figure 16 correspond to the mean parameters (b1, b2, b3) and deviation parameters (sd1, sd2, sd3), respectively. When the pseudorandom sequences were uncorrelated, mean bias was 0. With correlations of 0.1, 0.2, and 0.3 among pseudorandom sequences, bias in b1 reached 8%, 16%, and 24%, respectively. Larger correlations sometimes resulted in much greater bias. Naturally, MSE followed a similar pattern; 95% CI coverage started at 94% to 96% and decreased drastically with increasing correlation among pseudorandom sequences. Coverage for b1 was 77%, with a correlation of only 0.2. Once the correlation among pseudorandom sequences reached 0.7, bias, MSE, and coverage improved somewhat for mean estimates but continued to worsen for deviation estimates. With the correlation set at 0.9, 1 model out of 1000 failed to produce a standard error for sd1, and another failed to produce standard errors for b1 through b3 and sd1 through sd2.

Figure 16

Effects of Correlated Random Parameters on Bias, Absolute Bias, MSE, and 95% CI Coverage.

Method

To examine the joint effect of the number of random parameters and the number of Halton draws on parameter estimates and standard errors, we ran a series of RPL models with the actual (not simulated) data from a DCE on preferences for allocation of organs for transplant.³⁷ In the original study, each of 2051 respondents completed 30 choice tasks with 2 alternatives and 15 attributes. The analysis variables included 24 attribute-related parameters and 1 alternative-specific constant. Based on published health-related DCEs (Figure 1), we selected 3 levels for number of random parameters (5, 10, and 15) and 8 levels for number of draws (250, 500, 750, 1000, 2500, 5000, 7500, and 10 000). Because some results did not stabilize with 10 000 draws, we implemented 4 additional levels for number of draws (12 500, 15 000, 17 500, and 20 000). Using the R gmnl package, we ran RPL models (not generalized MNL, as the reader may infer from the package name) with all permutations of number of random parameters and number of draws, assuming normally distributed random parameters (Appendix B: Table B1). We arbitrarily selected as random effects the first 5, 10, or 15 attribute-related parameters. For each combination of number of random parameters and number of draws, we ran a single model.

As an example, the model below specified the first 5 attribute-related parameters as random:

where U_nk = utility for person n and alternative k; asc1 is an alternative-specific constant for alternative 1; n = 1, 2, …, N; β_pn is a normally distributed random parameter associated with attribute p; β_q is a population-level parameter associated with attribute q; and ε_nk is the independently and identically distributed random error, following the Gumbel distribution.

We compared bias and the precision of results across models. For each number of random parameters, we treated the model with the maximum number of draws (20 000) as the gold standard. To assess bias in estimated means, we calculated the MRS by dividing each parameter estimate by the parameter estimate for posttransplant quality of life (QOL). (For posttransplant QOL, we divided by the parameter estimate for pretransplant QOL.) We compared each MRS to the corresponding MRS from the maximum-draws model. To assess bias in estimated deviations, we calculated the coefficient of variation (COV)—the SD estimate divided by the corresponding mean estimate—and compared each COV to the corresponding COV from the maximum-draws model. To assess precision, we calculated relative efficiency by dividing the RSE (SE divided by parameter estimate) from the maximum-draws model by the corresponding RSEs from other models. A percentage <100% would indicate that the current estimate is less efficient than the gold-standard estimate because the current estimate has a higher RSE. We expressed bias and relative efficiency as percentages relative to the estimates from the maximum-draws model.

We also recorded the central processing unit (CPU) time, memory use, and log-likelihood (LL) for each model and plotted these values by number of random parameters and number of draws.

Results

For the real-data models, Table 2 summarizes the degree of bias (top half of table) and the relative efficiency (bottom half of table) compared with the estimates from the maximum-draws model by number of random parameters and parameter type across models with different numbers of draws. Models with 5 random parameters rarely showed substantial bias. Models with 10 random parameters sometimes exhibited bias, especially in estimated deviations. For models with 15 random parameters, bias often exceeded 50%, even for estimated means with no corresponding random effect, and the most biased estimate differed from the maximum-draws estimate by a factor of 362 (maximum, 36 238.7). Nineteen models, all with 15 random parameters, had at least 1 estimate with bias >1000% in absolute value (not shown). Ten of those 19 models used >1000 draws, and 2 used 17 500 draws, indicating that parameter estimates in approximately 15 random-parameter models had not yet stabilized by the time the number of draws reached 20 000.

Table 2

Absolute Value of Percent Bias and Relative Efficiency by Number of Random Parameters Estimated and Parameter Type.

With regard to relative efficiency (ratio of standard error to estimate in the maximum-draws model divided by the same ratio from the current model), in models with 5 random parameters, the RSEs for estimated means differed little from the RSEs in the maximum-draws model. The relative efficiency of estimated deviations exhibited some variation. Models with 10 random parameters produced a similar pattern of results but much greater variation in the relative efficiency of estimated deviations. With 15 random parameters, relative efficiency fluctuated greatly for all types of parameter estimates. For the models with 15 random parameters, 16% (54/330) of the relative efficiencies for random parameters and their deviations could not be estimated because the corresponding model did not produce a standard error estimate. The number of draws in these problematic models ranged from 250 to 20 000. The large relative efficiencies with smaller numbers of draws, mostly for deviation parameters, suggest that the stability of deviation parameters was especially sensitive to the number of draws used.

Figure 17 shows how the bias (panels a1, b1, and c1) and relative efficiency (panels a2, b2, and c2) of parameter estimates varied by number of draws (x-axis) across models with 5, 10, or 15 random parameters (panels a1-a2, b1-b2, and c1-c2, respectively; for visibility, the y-axis scale varies across plots.) Each point in the figure represents a single parameter estimate from a single model. Estimates became increasingly stable as the number of draws increased. In models with 5 random parameters, bias and the relative efficiency of estimated means (panels a1 and a2, plots 1-5) varied little but were more stable with ≥1000 draws. For certain estimated deviations, greater bias and efficiency differences emerged (panels a1 and a2, plot 6). Bias and the relative efficiency of estimated deviations stabilized after 1000 to 2500 draws.

Figure 17

Bias in Parameter Ratios and Relative Efficiency of Parameter Estimates by Number of Random Parameters Estimated and Number of Halton Draws.

In models with 10 random parameters, the bias and relative efficiency of most estimated means (Figure 17, panels b1 and b2, plots 1-5) varied little but appeared slightly more stable after about 1000 to 2500 draws. In 1 exceptional case (panels b1 and b2, plot 4), bias and relative efficiency did not stabilize before the number of draws reached 17 500 (bias near 5% and relative efficiency near 105%) and may not have stabilized at all. For some estimated deviations (panels b1 and b2, plots 6 and 7), estimates clearly did not stabilize even as the number of draws approached 20 000.

Models with 15 random parameters (Figure 17, panels c1 and c2) showed substantial bias, relative efficiency differences, and fluctuation in bias and relative efficiency. These results emerged for all types of parameters, including estimated means with no corresponding random effect. Results stabilized somewhat after 1000 draws but did not completely stabilize even with 20 000 draws.

The required amount of computer memory was a linear function of the number of Halton draws, with 12 to 13 GB required for 1000 draws, 27 to 29 GB required for 2500 draws, and 221 to 242 GB for 20 000 draws (not shown). CPU time was a quadratic function of the number of draws, with 6 to 7 hours required for 1000 draws, 9 to 11 hours for 2500 draws, and 92 to 122 hours (4-5 days) for 20 000 draws. LL appeared to stabilize after about 2500 draws with 5 random parameters and after 15 000 draws with 10 random parameters but never clearly stabilized with 15 random parameters (not shown).

Methods

Using R version 3.5.1, we simulated 500 study samples with n = 500 individuals per sample. The sample data included 3 normally distributed random parameters (b1, b2, b3) with means 1/3, 2/3, and 1.0 and SDs 1/9, 2/9, and 1/3 respectively. The data-generation model was ${\text{U}}_{\text{nk}}={\text{β}}_{\text{1n}}\text{*x}{1}_{\text{k}}+{\text{β}}_{2\text{n}}\text{*x}{2}_{\text{k}}+{\text{β}}_{3\text{n}}\text{*x}{3}_{\text{k}}+{\text{ε}}_{\text{nk}}$, where U_nk = utility for person n and alternative k; n = 1, 2, …, N; and ε_nk is the random error, following the Gumbel distribution.

We designed a simple DCE with 6 choice tasks, 3 alternatives per task, and 3 attributes (corresponding to b1, b2, and b3) per alternative. Then, for each of the 500 samples, we generated 2 sets of bootstrapped samples with replacement at 100% of the original sample size. One set of bootstrapped samples used 1000 replications, and the other used 2500 replications. We ran MNL models (Appendix B: Table B1) using the original sample and each set of bootstrapped samples as well as an RPL model using the original sample. For simplicity and ease of comparison, we generated the data and estimated the models without intercepts. To obtain bootstrapped estimates for each parameter, we calculated the mean of the parameter estimate across each set of bootstrapped samples. To estimate bootstrapped confidence intervals, we used 3 methods: (1) the percentile method (2.5th and 97.5th percentiles of the parameter estimate across bootstrapped samples); (2) the standard deviation method (SD of the parameter estimate across bootstrapped samples); and (3) the bias-corrected and accelerated method, which is similar to the percentile method but adjusts the percentiles to account for bias (the percentage of bootstrapped estimates that are less than the conventional estimate) and acceleration (sensitivity of the bootstrapped standard error to the true value of the parameter being estimated).⁴³ We then compared the estimates to the true parameter values with respect to mean percent bias (the average difference between the estimate and the true parameter value, expressed as a percentage of the true parameter value), mean absolute percent bias (the average absolute difference between the estimate and the true parameter value, expressed as a percentage of the true parameter value), MSE (the mean of the square of bias), and 95% CI coverage (the percentage of samples in which the 95% CI included the true parameter value). We expected that, compared with the conventional MNL model, bootstrapping would yield confidence interval coverage closer to the nominal (95%) coverage.

As described previously in the sections on stakeholder engagement and changes to the research protocol, these methods reflect simplification of the initial study protocol to allow completion of additional study components. Specifically, we separated the bootstrapping analysis from the simulations exploring the effects of DCE design features and model assumptions on DCE results. In addition to reducing the computing resources required for the study, this change enabled us to conduct a thorough comparison of bootstrapped vs traditional methods without confounding the effects of the variance estimation method with those of DCE design features or model assumptions.

Results

Table 3 shows the results of the bootstrapping simulations. The bootstrapped parameter estimates had little bias (mean percent bias, −4.5% to 0.6%). The MSE for the bootstrapped estimates of b1 and b2 was similar to the MSE for the RPL estimates, but the MSE for the bootstrapped estimate of b3 was higher than the MSE for the RPL estimate. Confidence interval coverage for b1 was close to the nominal value of 95% regardless of estimation method, but only the RPL model provided adequate coverage for all 3 parameters, and coverage for b3 was <75% for the conventional and bootstrapped MNL models. In fact, as the true parameter means and SDs increased from b1 to b2 to b3, bias in the MNL estimates increased, and the confidence interval coverage of those estimates (regardless of bootstrapping method) decreased.

Table 3

Results of Simulations Assessing Bootstrapped MNL Estimation of Parameter Means in the Presence of Preference Heterogeneity.

Substudy Limitations

One limitation of this component of the study is its complexity. We were able to conduct simulations based on only 2 actual DCEs, and we reduced the number of simulation parameters as well as the number of levels for certain parameters in the parent DCEs. Nonetheless, we obtained a massive quantity of simulated data. We compared selected simulation scenarios with a limited number of default scenarios (pilot sample size of 240 participants, no unmeasured interaction [if all attributes were included in the data-generating process], variance decreased by 10% [with selected attributes excluded], no selection bias, no ignoring of attributes, main DCE sample size of 1000 participants, full analysis model or main DCE, unmeasured interaction in the data-generating process when comparing across different main DCE analysis models). Although we believe that the comparisons and default scenarios were well chosen, we may have merely scratched the surface of the available information.

Related to the complexity of this study component, our aim was descriptive rather than explanatory, which limits the degree to which we can describe the mechanisms underlying the observed patterns and discuss the practical implications. Another limitation is that our findings and interpretations may apply only to the specific simulation scenarios we examined. For example, the data-generating mechanism included only normally distributed random parameters, and some of our simulation parameters—especially those not as strongly supported by research evidence (eg, the amount of selection bias, 50% of the population ignoring certain attributes)—may have been unrealistic. A related limitation involves the selection of real DCEs on which to base the simulations. Using data from studies typical of recent health-related DCEs would be ideal. However, we believed that replication of the original results was an important step in re-creating the original analysis models, and we were able to access and transmit data from only a few studies. Although study 2 was typical except for the large number of choice tasks, study 1 was atypical (ie, it had a large sample, unusual study focus, and many choice tasks and attributes). To the extent that these studies fail to represent recent health-related DCEs, our conclusions may have limited applicability. Also, any flaws in the original studies may have affected the sets of parameter estimates on which we based our simulations. On the other hand, our simulations were much more realistic than they would have been if we had simply invented population parameters.

Summary

Researchers may take comfort from the fact that in our simulation settings, pilot study problems had no general effect on main DCE results. However, our findings also suggest that unusual DCE designs (eg, those with highly correlated attributes) may require an abundance of caution (eg, thorough sensitivity analyses). Research should continue to explore the compound effects of pilot and main DCE design features and assumption violations on findings by searching both more deeply (eg, those with more analyses) and more broadly (eg, those with different sets of simulation parameters).

Limitations

This component of our study has at least 2 important limitations. First, coverage of the 0 to 1 interval, correlations among random parameters, and significance of normality tests do not translate directly to the validity and reliability of findings. However, we used simulation to demonstrate the effects of correlated random parameters on bias and confidence interval coverage, and we used real data to examine how the number of random parameters and number of draws affect the validity and reliability of estimates. Second, our real-data example came from a single study, and our simulations focused on specific parameter combinations; our findings and interpretations are limited to the specific data set and simulation scenarios we examined. Replication with different data sources and simulation parameters would provide more context for our findings, interpretations, and conclusions. Specifically, although we suspect that health-related DCEs should use far more Halton draws than they have in the past and we know that adding random parameters requires increasing the number of draws, our findings do not provide precise thresholds by which to judge DCEs.

Summary

Despite these limitations, given our findings and the current state of health-related DCEs, we make the following recommendations:

• The number of random parameters being estimated should inform the number of Halton draws to use.
• RPL analysis should include a sensitivity analysis to assess the stability of estimates with larger numbers of draws, especially when estimating many random parameters.
• Reports of RPL analyses should include the number of random parameters estimated, the number of draws used, and the sensitivity analysis results.
• Researchers should develop more detailed guidelines for RPL analyses.
• Researchers should continue to explore alternatives to Halton draws for simulated maximum likelihood estimation.

With regard to the last recommendation, randomized and shuffled Halton sequences do not solve the problems we have observed, but Sobol sequences may improve on the performance of Halton sequences, and modified Latin hypercube sampling—the current standard in the transport field—also deserves more attention.⁷ Pseudorandom sampling might prevent the problems we observed but would greatly increase the amount of computing time required.²³

For the analyst who chooses to use Halton draws, because our simulation results are specific to our simulation settings and parameters, we have recommended no specific numeric thresholds for number of random parameters or number of draws. These thresholds should be decided by the analyst. However, our analyses describing violations of the independence and normality assumptions, which were not simulations but rather examinations of Halton sequences, clearly indicated that estimating more than 5 to 7 normally distributed random parameters without using thousands of draws would lead to such violations. In our simulations examining bias and confidence interval coverage with 3 correlated, normally distributed random parameters, we found that in our setting, minor assumption violations in the form of correlations of 0.1 or 0.2 led to substantial estimation problems. One analytic approach would be to decide on acceptable thresholds by combining these findings, run initial analyses using those thresholds, and follow up with a sensitivity analysis using the methods from our real-data example (essentially convergence plots of estimates and SEs for different numbers of draws). The set of acceptable thresholds should incorporate the knowledge that the number of draws should increase more quickly than the square root of the number of observations (individuals × choice tasks).^7,23 This strategy may provide a realistic way to manage the observed resource constraints without sacrificing study quality. The strategy could be enhanced by using software that provides warnings regarding the number of draws used⁴⁸ or allows comparison of multiple methods.⁴⁹

Table A1. Study 1 attributes and variables (PDF, 36K)

Table A2. Study 2 attributes and variables (PDF, 34K)

Figure A1. Study 1 convergence plots, “best case” scenario (PDF, 714K)

Figure A2. Study 2 convergence plots, “best case” scenario (PDF, 498K)

Figure A3. Study 1 convergence plots, “worst case” scenario (PDF, 647K)

Figure A4. Study 2 convergence plots, “worst case” scenario (PDF, 452K)

Figure A5. Percent bias in ratios of parameter estimates by pilot simulation parameters and main DCE sample size, Study 1 (PDF, 757K)

Figure A6. Percent bias in ratios of parameter estimates by pilot simulation parameters and main DCE sample size, Study 2 (PDF, 563K)

Figure A7. Percent bias in ratios of parameter estimates by pilot simulation parameters and main DCE analysis model, Study 1 (PDF, 783K)

Figure A8. Percent bias in ratios of parameter estimates by pilot simulation parameters and main DCE analysis model, Study 2 (PDF, 526K)

Figure A9. Relative standard error by pilot simulation parameters and main DCE sample size, Study 1 (PDF, 928K)

Figure A10. Relative standard error by pilot simulation parameters and main DCE sample size, Study 2 (PDF, 525K)

Figure A11. Relative standard error by pilot simulation parameters and main DCE analysis model, Study 1 (PDF, 716K)

Figure A12. Relative standard error by pilot simulation parameters and main DCE analysis model, Study 2 (PDF, 496K)

Figure A13. D-error by pilot simulation parameters and main DCE sample size, Study 1 (PDF, 123K)

Figure A14. D-error by pilot simulation parameters and main DCE sample size, Study 2 (PDF, 139K)

Figure A15. D-error by pilot simulation parameters and main DCE analysis model, Study 1 (PDF, 114K)

Figure A16. D-error by pilot simulation parameters and main DCE analysis model, Study 2 (PDF, 121K)

Table B1. Analysis models (PDF, 147K)