Flexibility when generating Bayesian efficient designs

[Bob123]

Dear Michiel,

When generating an efficient design, I understand that the process looks roughly something like: starting position e.g., a D-Efficient design, conduct pilot in ~10% of your sample using this design, run the model (e.g., an MNL model in Apollo) using the same model spec as specified in ngene, and finally use betas and SE's as priors to generate the Bayesian efficient design.

However, for best practice (or in real-practice) how flexible is the point between obtaining priors and generating your second (the Bayesian efficient) design and how often does the specification remain exactly the same or change? Is it ok for example, to remove an interaction in the efficient design because you saw in initial modelling that is was unlikely to be significant. Also, what is the protocol for if something unexpected e.g., a beta of an important main effect being (-ve) rather than (+ve) (and definitely wrong) occurs? How often would you "tweak" priors, or only use some of them in your Bayesian efficient design?

Specifically, I have an unlabelled experiment with lots of categorical scenario/ contextual variables that I included in my initial design on 10% of my sample as interactions. Once I ran the model however, one of my main effects is the wrong sign (+ instead of -) but corrects when I remove some of the interactions? For my next step i.e., adding priors into ngene and going from d-efficient to bayesian efficient, is it ok to remove some interactions and maybe even fiddle with the signs of some coefficients or their magnitude? For reference, my pilot was n=70.

In your course I remember you saying "things can go badly wrong if you misspecify your priors" which I why I would like to remove some that are non-sensical and seem to make others also make more sense. I also remember you saying that guessing priors or fiddling around should be left to the experts!

Thanks in advance for any advice.

Rob

[Michiel Bliemer]

You have a lot of flexibility because it is all about best guesses for utility function and priors. So as long as your choices are not unreasonable it is fine.

You can take out an interaction effect (so you estimate a model in Apollo without the interaction effect to get priors).

If a parameter has the wrong sign, then I would manually correct it to be the correct sign, for example a uniform distribution bounded from above or below by 0.

Things can go really wrong if your priors are too far from 0, so if you would use b1[0.5] * X[100,200,300], where 0.5 is far too large for the large numerical levels. So as long as you do not choose some outrageously large priors then you should be ok. In case of categorical variables, where you multiply with 0 or 1 with dummy coding you would typically expect parameter values to be in the range of [-1,1].

Michiel

[Bob123]

Dear Michiel,

Thanks for you prompt reply.

As you suggested, I removed some interactions that I had in my original ngene model specification and ran the model (in Apollo) to obtain priors. Below is the output:

Code: Select all

Estimates:
                                       Estimate        s.e.   t.rat.(0)  p(1-sided)    Rob.s.e. Rob.t.rat.(0)  p(1-sided)
asc_DrugA                              0.000000          NA          NA          NA          NA            NA          NA
b_CERT_very_low                        0.000000          NA          NA          NA          NA            NA          NA
b_CERT_low                             0.404195    0.276399      1.4624     0.07182    0.296966        1.3611    0.086744
b_CERT_moderate                        1.206697    0.320756      3.7620   8.427e-05    0.372031        3.2435  5.9028e-04
b_CERT_high                            1.705627    0.450533      3.7858   7.661e-05    0.554702        3.0749    0.001053
b_WAIT                                -0.028868    0.041982     -0.6876     0.24584    0.042642       -0.6770    0.249211
b_LIFExWAIT                          1.5080e-04    0.001103      0.1368     0.44561    0.001071        0.1407    0.444037
b_QOL_light_workxWAIT                  0.011200    0.012144      0.9222     0.17820    0.011151        1.0044    0.157600 [?]
b_QOL_comp_disabledxWAIT               0.014339    0.012937      1.1084     0.13385    0.012188        1.1765    0.119696 [?]
b_BENEFIT_smallxWAIT                  -0.024892    0.013533     -1.8394     0.03293    0.015589       -1.5967    0.055161 [?]
b_BENEFIT_substantialxWAIT             0.002079    0.011952      0.1740     0.43094    0.012302        0.1690    0.432881 [?]
b_QOL_light_workxCERT_very_low        -0.096430    0.356192     -0.2707     0.39330    0.315712       -0.3054    0.380017 [?]
b_QOL_comp_disabledxCERT_very_low      0.447222    0.340265      1.3143     0.09437    0.325787        1.3727    0.084916 [?]
b_QOL_light_workxCERT_high             0.406101    0.348976      1.1637     0.12227    0.351122        1.1566    0.123722 [?]
b_QOL_comp_disabledxCERT_high         -0.248969    0.316094     -0.7876     0.21545    0.374766       -0.6643    0.253239 [?]
b_BENEFIT_smallxCERT_very_low         -0.378522    0.382386     -0.9899     0.16111    0.403617       -0.9378    0.174167 [?]
b_BENEFIT_substantialxCERT_very_low   -0.383123    0.354293     -1.0814     0.13977    0.389200       -0.9844    0.162463 [?]
b_BENEFIT_smallxCERT_high             -0.117946    0.348369     -0.3386     0.36747    0.383841       -0.3073    0.379316 [?]
b_BENEFIT_substantialxCERT_high       -0.767797    0.346004     -2.2190     0.01324    0.335874       -2.2860    0.011128 [?]


In the output above, I am quite happy with my beta's for main effects (they are the expected signs (+ve/-ve) and magnitude relative to each other (preference order makes sense). I have included these in my ngene Bayesian efficient design as normally distributed Bayesian priors (b1 below). I am also mostly happy with the estimates for one of the interactions (LIFExWAIT) as it is the expected sign, as you have suggested in previous posts, I have included this as a uniform distributions Bayesian prior with a lower bound of 0 so ngene doesn't take draws with a negative estimate (b2 below).

With the remaining interactions (marked with a [?] in my output above), I am not sure what is the best way to proceed. Because they are interactions between categorical variables with multiple levels, it is difficult to interpret whether they should be +ve/-ve and I don't think it would be appropriate to at this stage because they are not straightforward.

*To avoid making an error in adding Bayesian priors, would you agree that the following approach would work for those which I am uncertain of:*
- For those interactions with small beta's close to 0 (+ve or -ve), input as normally distributed priors around 0 with a quite large SE estimate e.g. 0.3 so ngene considers a wide variety of draws for these parameters? - see i6-i8 below.
- For the interactions with quite large beta's (e.g., with BENEFIT - the bottom four (i12 to i15 below) are all negative between -0.3 and -0.7) even though I don't know whether these estimates are reliable, include them anyway? They have not changed much in modelling when removing interactions (always same sign and quite large in magnitude) so hope for the best? Is there a way I could make this approach a bit less risky?

For reference, I have tried out my suggested approach in ngene and is coded as follows. I seem to get reasonable sp estimates and can't see anything immediately wrong with the choice tasks it creates

Code: Select all

 
;model(m2):
U(DrugA) = c1
         + b1.dummy[(n,0.4,0.3)|(n,1.2,0.4)|(n,1.7,0.6)]   *  CERT[1,2,3,0]
         + b2[n,-0.03,0.04]                                *  WAIT[0,0.5,1,2]
 
         + i3[(u,-0.005,0)]                               *  WAIT * LIFE[0.5,1,2,3]
         + i4[(n,0.01,0.01)]                              *  WAIT * QOL.dummy[0] 
         + i5[(n,0.01,0.01)]                              *  WAIT * QOL.dummy[3] 
         + i6[(n,0,0.3)]                                  *  WAIT * BENEFIT.dummy[0]
         + i7[(n,0,0.3)]                                  *  WAIT * BENEFIT.dummy[2]
         + i8[(n,0,0.3)]                                  *  CERT.dummy[0] * QOL.dummy[0] 
         + i9[(n,0.45,0.33)]                              *  CERT.dummy[0] * QOL.dummy[3] 
         + i10[(n,0.4,0.35)]                              *  CERT.dummy[3] * QOL.dummy[0] 
         + i11[(n,-0.25,0.37)]                            *  CERT.dummy[3] * QOL.dummy[3]
         + i12[(n,-0.38,0.4)]                             *  CERT.dummy[0] * BENEFIT.dummy[0]
         + i13[(n,-0.38,0.39)]                            *  CERT.dummy[0] * BENEFIT.dummy[2]
         + i14[(n,-0.12,0.38)]                            *  CERT.dummy[3] * BENEFIT.dummy[0]
         + i15[(n,-0.77,0.34)]                            *  CERT.dummy[3] * BENEFIT.dummy[2]


Thanks in advance again!

Rob


Full ngene code with restrictions etc in case useful:

Code: Select all

design
;alts(m1) = DrugA*, DrugB*
;alts(m2) = DrugA*, DrugB*
;rows = 40
;block= 4
;eff  = m2(mnl,d,median)
;bdraws = halton(200)
;alg = mfederov
;reject:
         DrugA.LIFE = 0.5 and DrugA.WAIT = 0.5,
         DrugA.LIFE = 0.5 and DrugA.WAIT = 1,
         DrugA.LIFE = 0.5 and DrugA.WAIT = 2,
         DrugA.LIFE = 1   and DrugA.WAIT = 1,
         DrugA.LIFE = 1   and DrugA.WAIT = 2,
         DrugA.LIFE = 2   and DrugA.WAIT = 2,
         DrugB.LIFE = 0.5 and DrugB.WAIT = 0.5,
         DrugB.LIFE = 0.5 and DrugB.WAIT = 1,
         DrugB.LIFE = 0.5 and DrugB.WAIT = 2,
         DrugB.LIFE = 1   and DrugB.WAIT = 1,
         DrugB.LIFE = 1   and DrugB.WAIT = 2,
         DrugB.LIFE = 2   and DrugB.WAIT = 2
;require:
         DrugA.LIFE = DrugB.LIFE,
         DrugA.QOL  = DrugB.QOL,
         DrugA.BENEFIT  = DrugB.BENEFIT
;model(m1):
U(DrugA) = c1
         + b1.dummy[0|0|0]                  *  QOL[1,2,3,0]
         + b2                               *  LIFE[0.5,1,2,3]
         + b3.dummy[0.0001|0.0002]          *  BENEFIT[1,2,0]
         + b4.dummy[0.0001|0.0002|0.0003]   *  CERT[1,2,3,0]
         + b5[-0.0001]                      *  WAIT[0,0.5,1,2]
 
         /
U(DrugB) = b1.dummy[0|0|0]                  *  QOL[1,2,3,0]
         + b2                               *  LIFE[0.5,1,2,3]
         + b3.dummy[0.0001|0.0002]          *  BENEFIT[1,2,0]
         + b4.dummy[0.0001|0.0002|0.0003]   *  CERT[1,2,3,0]
         + b5[-0.0001]                      *  WAIT[0,0.5,1,2]
;model(m2):
U(DrugA) = c1
         + b1.dummy[(n,0.4,0.3)|(n,1.2,0.4)|(n,1.7,0.6)]   *  CERT[1,2,3,0]
         + b2[n,-0.03,0.04]                                *  WAIT[0,0.5,1,2]
 
         + i3[(u,-0.005,0)]                               *  WAIT * LIFE[0.5,1,2,3]
         + i4[(n,0.01,0.01)]                              *  WAIT * QOL.dummy[0] 
         + i5[(n,0.01,0.01)]                              *  WAIT * QOL.dummy[3] 
         + i6[(n,0,0.3)]                                  *  WAIT * BENEFIT.dummy[0]
         + i7[(n,0,0.3)]                                  *  WAIT * BENEFIT.dummy[2]
         + i8[(n,0,0.3)]                                  *  CERT.dummy[0] * QOL.dummy[0] 
         + i9[(n,0.45,0.33)]                              *  CERT.dummy[0] * QOL.dummy[3] 
         + i10[(n,0.4,0.35)]                              *  CERT.dummy[3] * QOL.dummy[0] 
         + i11[(n,-0.25,0.37)]                            *  CERT.dummy[3] * QOL.dummy[3]
         + i12[(n,-0.38,0.4)]                             *  CERT.dummy[0] * BENEFIT.dummy[0]
         + i13[(n,-0.38,0.39)]                            *  CERT.dummy[0] * BENEFIT.dummy[2]
         + i14[(n,-0.12,0.38)]                            *  CERT.dummy[3] * BENEFIT.dummy[0]
         + i15[(n,-0.77,0.34)]                            *  CERT.dummy[3] * BENEFIT.dummy[2]
 
         /
U(DrugB) = b1.dummy[(n,0.4,0.3)|(n,1.2,0.4)|(n,1.7,0.6)]   * CERT[1,2,3,0]
         + b2[n,-0.03,0.04]                                *  WAIT[0,0.5,1,2]
 
 
         + i3[(u,-0.005,0)]                               *  WAIT * LIFE[0.5,1,2,3]
         + i4[(n,0.01,0.01)]                              *  WAIT * QOL.dummy[0] 
         + i5[(n,0.01,0.01)]                              *  WAIT * QOL.dummy[3] 
         + i6[(n,0,0.3)]                                  *  WAIT * BENEFIT.dummy[0]
         + i7[(n,0,0.3)]                                  *  WAIT * BENEFIT.dummy[2]
         + i8[(n,0,0.3)]                                  *  CERT.dummy[0] * QOL.dummy[0] 
         + i9[(n,0.45,0.33)]                              *  CERT.dummy[0] * QOL.dummy[3] 
         + i10[(n,0.4,0.35)]                              *  CERT.dummy[3] * QOL.dummy[0] 
         + i11[(n,-0.25,0.37)]                            *  CERT.dummy[3] * QOL.dummy[3]
         + i12[(n,-0.38,0.4)]                             *  CERT.dummy[0] * BENEFIT.dummy[0]
         + i13[(n,-0.38,0.39)]                            *  CERT.dummy[0] * BENEFIT.dummy[2]
         + i14[(n,-0.12,0.38)]                            *  CERT.dummy[3] * BENEFIT.dummy[0]
         + i15[(n,-0.77,0.34)]                            *  CERT.dummy[3] * BENEFIT.dummy[2]
$
 
 

[Bob123]

*EDIT*

Apologies for an extra question.

In the above post I mentioned that I can't see anything immediately wrong with the choice tasks the ngene code below creates.

I realise I use quite a few restrictions so as to not present dominant alternatives and also keep attributes constant (LIFE, QOL and BENEFIT) across alternatives. But in the choice sets generated by the code (in any of the sets produced), there are no 0.5 levels for LIFE which I would expect to be presented at least a few times when WAIT is 0 - the restrictions on the model should not prevent this? Also in 40 choice tasks, LIFE = 2 only comes up 3 times. Perhaps this, (in addition to a smaller sample size for pilot study) is why some of my betas are not as impactful on utility as predicted.

In the code, could you possibly identify why LIFE=0.5 is not occurring in the choice sets and LIFE=2 only occurs 3 times, should this be the case or should I change something so these combinations are presented to respondents?

Thanks again so much,

Rob

[Michiel Bliemer]

Hi Rob,

These questions are going a bit beyond the support I can offer as these are very specific questions about a particular study. But to give some guidance:

- Your way of setting priors looks reasonable.

- Your b2 is now a random parameter, you need to make it a Bayesian prior with brackets around the prior.

- When inspecting a design, you can look in the output at Factorial, which shows you the choice tasks that Ngene uses in the candidate set for the modified Federov algorithm. Here you can see that LIFE = 0.5 never appears, so this means that your constraints do not allow this level. This may have to do with some of your constraints, or with the dominance constraints that Ngene applies. You could remove the dominance checks in Ngene by removing the asterisk (*) after the alternative names, or you should relax the reject constraints that you impose.

- To ensure that you a sufficient number of levels, you could impose attribute level constraints such as using LIFE[0.5,1,2,3](0-40,0-40,5-40,0-40) to ensure that level 2 appears at least 5 times. Note that each constraint makes it more difficult for Ngene to generate designs that fit all constraints and you will lose some efficiency.

- You have a very large number of Bayesian priors so you need different types of draws and much more of them. I suggest using ;bdraws = sobol(1000)

- Because you have so many constraints, only 1440 choice tasks satisfy your constraints. It is a very computationally demanding script, so I suggest using ;alg = mfederov(candidates = 1000)

Michiel

[Bob123]

Dear Michiel,

Thanks so much for taking the time to answer these questions - of course, your expert advice worked great and I was able to generate my final Bayesian efficient design!

Best,

Rob