How to incorporate weights into the likelihood of a WinBUGS model?

Question

I want to incorporate weights into the likelihood of a WINBUGS model to do what brms does with weights.

The usual BUGS approaches to accomplish that for dnorm and dpois are not working for dbin.

In brms, as paul.buerkner says here, that is accomplished with a Stan code like this:

vector[N] weights; \\ model weights 
target += weights[n] * neg_binommial_2_log_lpmf(Y[n] | mu[n], shape);

When I implement this approach in my BUGS model I get the error below detailed.

Here are the data and model:

library(R2WinBUGS)
dat <- data.frame(
  A = c(1, 1, 0, 0), B = c(1, 0, 1, 0),
  Pass = c(278, 100, 153, 79), Fail = c(743, 581, 1232, 1731), Weights= c(3, 1, 12, 3))
N <- length(dat$Pass)
case <- dat$Pass
nn <- dat$Fail+dat$Pass
A <- dat$A
B <- dat$B
data <- list (N = N, case = case, A = A, B = B, nn = nn)
weights <- dat$Weights
data.weighted <- list (N = N, case = case, A = A, B = B, nn = nn, weights = weights)
model {
for (i in 1:N){
pp[i] <- exp(a0) * (1 + a1*A[i] + a2*B[i] + a3*A[i]*B[i])
case[i] ~ dbin(pp[i], nn[i])}
a0 ~ dnorm(0, 5)
a1 ~ dnorm(0, 5)
a2 ~ dnorm(0, 5)
a3 ~ dnorm(0, 5)
ones <- 1
ones ~ dbern(C1)
C1 <- step(0 - a0)*step(1 + a1)*step(exp(- a0) - 1 - a1)*step(1 - a2)*
step(exp(-a0) -1 - a2)*step(1 + a1 + a2 + a3)*step(exp(- a0) - 1 - a1 - a2 - a3)}
inits <- function(){list(a0=-4,a1=0,a2=0, a3=0)}
sim1 <- bugs(data, inits=inits, model.file = "C:/.../bug/regi1.bug",
parameters = c("a0", "a1", "a2", "a3"),
n.chains = 4, n.burnin = 500000, n.iter = 1000000,
bugs.directory = "C:/.../winbugs14_unrestricted/WinBUGS14", debug=FALSE)
print(sim1)

When I implement the solution like that used in brms, I get the error:

Incompatible copy

 BugsCmds.TextError   [000003A1H]   .beg    INTEGER 1670313368
    .end    INTEGER 44779932    .name   ARRAY 256 OF
CHAR    "C:/Users/.../AppData/Local/"   ...     .pos    INTEGER 144
    .text   TextModels.Model    [02AB7940H]     .v  Views.View  [02AB8C80H]  
BugsCmds.Parse   [00000470H]    .name   ARRAY 256 OF
CHAR    "C:/Users/.../AppData/Local/"   ...
    .s  BugsScanners.Scanner    Fields  StdInterpreter.CallProc   [0000047AH]
    .a  BOOLEAN FALSE   .b  BOOLEAN FALSE   .c  BOOLEAN FALSE
    .i  Meta.Item   Fields  .imported   ARRAY 256 OF CHAR   ""   ...
    .importing  ARRAY 256 OF CHAR   ""   ...    .mn Meta.Name   "BugsCmds"
    .mod    StdInterpreter.Ident    "BugsCmds"  .object ARRAY 256 OF CHAR   ""  
...     .ok BOOLEAN TRUE    .parType    INTEGER 3   .pn Meta.Name   "Parse"
    .proc   StdInterpreter.Ident    "Parse"   ...   .res    INTEGER 0
    .v  StdInterpreter.ProcVal  Fields  .vi StdInterpreter.ProcIVal Fields
    .vii    StdInterpreter.ProcIIVal    Fields
    .vr StdInterpreter.ProcRVal Fields
    .vri    StdInterpreter.ProcRIVal    Fields
    .vrii   StdInterpreter.ProcRIIVal   Fields
    .vrr    StdInterpreter.ProcRRVal    Fields
    .vrri   StdInterpreter.ProcRRIVal   Fields
    .vrrii  StdInterpreter.ProcRRIIVal  Fields
    .vrs    StdInterpreter.ProcRSVal    Fields
    .vrsi   StdInterpreter.ProcRSIVal   Fields
    .vrsii  StdInterpreter.ProcRSIIVal  Fields
    .vs StdInterpreter.ProcSVal Fields
    .vsi    StdInterpreter.ProcSIVal    Fields
    .vsii   StdInterpreter.ProcSIIVal   Fields
    .vsr    StdInterpreter.ProcSRVal    Fields
    .vsri   StdInterpreter.ProcSRIVal   Fields
    .vsrii  StdInterpreter.ProcSRIIVal  Fields
    .vss    StdInterpreter.ProcSSVal    Fields
    .vssi   StdInterpreter.ProcSSIVal   Fields
    .vssii  StdInterpreter.ProcSSIIVal  Fields  StdInterpreter.Command  
[0000131CH]     .left   StdInterpreter.Ident    "BugsCmds"  .ptype  INTEGER 3
    .right  StdInterpreter.Ident    "Parse"   ... 
StdInterpreter.CallHook.Call   [00001441H]      .ch CHAR    0X  .e  ARRAY 64
OF CHAR ""   ...    .errorMsg   ARRAY 1 OF CHAR ""  .f  ARRAY 64 OF CHAR    ""
...     .g  ARRAY 64 OF CHAR    ""   ...
    .hook   StdInterpreter.CallHook [02B00050H]     .i  INTEGER 83
    .i0 INTEGER 0   .i1 INTEGER 0   .id StdInterpreter.Ident    "Parse"   ...
    .par0   Dialog.String   ""   ...    .par1   Views.Title ""   ...    .proc   ARRAY
240 OF CHAR "BugsCmds.Parse('C:/Users/..."   ...    .res    INTEGER 0
    .s0 Dialog.String   "C:/Users/.../AppData/Local/"   ...
    .s1 Dialog.String   ""   ...    .type   INTEGER 3   .x  INTEGER 0 
Dialog.Call   [00002FC8H]   .errorMsg   ARRAY 1 OF CHAR ""  .proc   ARRAY
240 OF CHAR "BugsCmds.Parse('C:/Users/..."   ...    .res    INTEGER 0 
BugsScript.Call   [00000130H]   .bugsCommands   ARRAY 240 OF
CHAR    "BugsCmds.SetOK; BugsCmds.Clear; "   ...    .i  INTEGER 114
    .item   Meta.Item   Fields  .j  INTEGER 82  .ok BOOLEAN FALSE
    .par    Dialog.Par  Fields  .pos    INTEGER -1  .res    INTEGER 0   .s  ARRAY 240
OF CHAR "BugsCmds.Parse('C:/Users/..."   ...    .scriptCommand  ARRAY 240
OF CHAR "#Bugs:check"   ...     .start  INTEGER 77
    .v  BugsScript.RECORD   Fields  BugsScript.Action.Do   [0000062FH] 
    .a  BugsScript.Action   [02C72400H]     .argNum INTEGER 0
    .bugsCommands   ARRAY 240 OF CHAR   "BugsCmds.SetOK; BugsCmds.Clear; "  
...     .p  ARRAY 3, 120 OF CHAR    Elements    .s  BugsScanners.Scanner    Fields
    .scriptCommand  ARRAY 240 OF CHAR   "#Bugs:check"   ...
    .vectorName BOOLEAN FALSE  Services.Exec   [00000136H] 
    .a  Services.Action [02C72400H]     .t  POINTER [67400170H] 
Services.IterateOverActions   [000002F4H] 
    .p  Services.Action [02C72400H]     .t  POINTER NIL
    .time   LONGINT 36194437  Services.StdHook.Step   [0000034DH] 
    .h  Services.StdHook    [02ABE380H]   HostWindows.Idle   [00004A86H] 
    .focus  BOOLEAN FALSE   .tick   Controllers.TickMsg Fields
    .w  HostWindows.Window  NIL  HostMenus.TimerTick   [00003422H] 
    .lParam INTEGER 0   .ops    Controllers.PollOpsMsg  Fields
    .wParam INTEGER 1   .wnd    INTEGER 4395250  Kernel.Try   [00003A61H] 
    .a  INTEGER 4395250     .b  INTEGER 1   .c  INTEGER 0
    .h  PROCEDURE   HostMenus.TimerTick  HostMenus.ApplWinHandler  
[00003841H]     .Proc   PROCEDURE   NIL     .hit    BOOLEAN FALSE
    .lParam INTEGER 0   .message    INTEGER 275     .res    INTEGER 0   .s  ARRAY 256
OF SHORTCHAR    "øëv"   ...     .w  INTEGER 1825694848  .wParam INTEGER 1
    .wnd    INTEGER 4395250 <system  (pc=7697BF1AH,  fp=0061FAB0H)
<system  (pc=769783E9H,  fp=0061FB98H) <system  (pc=76977C9DH, 
fp=0061FC14H) <system  (pc=7696CDCFH,  fp=0061FC1CH)  HostMenus.Loop
[00003BDEH]     .done   BOOLEAN FALSE   .f  SET {0..5}  .n  INTEGER 0
    .res    INTEGER 0   .w  HostWindows.Window  NIL  Kernel.Start   [00002B8CH]
    .code   PROCEDURE   HostMenus.Loop

Morris

I cannot test it myself but perhaps case[i] ~ weights[i] * dbin(pp[i], nn[i])} is all you need. – Riband

Thanks, @JuliusVainora. When I implement that solution, I get the error above (see EDIT above). Any thoughts? – Morris

I guess it means that the syntax is way off. No thoughts left, sorry. – Riband

Thanks, @JuliusVainora. – Morris

Related: stats.stackexchange.com/q/66832/10172. While in dnorm and dpois cases there is an intuitive way to introduce weights, it's somewhat more difficult with dbin. You can always add extra data according to your weights but, as the top answer there mentions, it's not a perfect solution. – Riband

Thanks, @JuliusVainora. Agree, add extra data is not the right way to go. Probably transferring the question to stats.stackexchange.com? Probably stats.stackexchange.com would be the right forum for this question. Any thoughts? – Morris

I guess adding extra data may be fine, but that depends on your ultimate goal. I'm not sure if it's possible to transfer the question while the bounty is still there. It would seem to be on-topic in both places, although clearly you are having no success here, so stats.se certainly makes sense. – Riband

I can vote to close the question and refer it to stats.se but (if it will work) another 4 votes will be needed then. Perhaps you could also flag your question and ask for a moderator to do that (since voting may take time). – Riband

Thanks, @JuliusVainora. I have flagged asking the moderator to transfer this question to stats.stackexchange.com. While it involves programming, the fact the usual BUGS tricks stats.stackexchange.com/questions/66832/… that work for dnorm and dpois are not working for dbin shows that transferring the question to a more statistical forum is justified. Please go ahead. – Morris

If the original dataset was collapsed to the unique rows, except that the weights variable indicates the number of times that row was observed, then multiplying a weight by a likelihood is an efficient way to calculate the likelihood of the parameters with respect to the original dataset. If, however, (as the OP suggests) the weights are some guess as to the number of unobserved people in the population with the same characteristics as a row observed in a sample, then weighting the likelihood this way is inconsistent with Bayes Rule because you are conditioning on unobserved outcomes. – Outstation

Thanks, @BenGoodrich. The weights variable in this [my] case aims to "[indicate] the number of times that row was observed, then multiplying a weight by a likelihood is an efficient way to calculate the likelihood of the parameters with respect to the original dataset." Just edited the question to make this clearer. BTW, any thoughts on how to incorporate the weights to the likelihood? – Morris

In Stan, you just multiply the weight value by the log likelihood. BUGS does not evaluate the likelihood so I think you have to loop over the weights. – Outstation

Thanks, @BenGoodrich. Could you kindly help me apply the "loop over the weights" in my BUGS code? Thank you in advance. – Morris

for (w in 1:weights[i]) case[i] ~ dbin(pp[i], nn[i]) – Outstation

Thanks @BenGoodrich. I hope I am getting it right but when I do model { for (i in 1:N){ pp[i] <- exp(a0) * (1 + a1*A[i] + a2*B[i] + a3*A[i]*B[i]) for (w in 1:weights[i]) { case[i] ~ dbin(pp[i], nn[i])} } a0 ~ dnorm(0, 5) a1 ~ dnorm(0, 5) a2 ~ dnorm(0, 5) a3 ~ dnorm(0, 5) # ensuring appropriate constraints for a0, a1, a2 and a3 ones <- 1 ones ~ dbern(C1) C1 <- step(0 - a0)*step(1 + a1)*step(exp(- a0) - 1 - a1)*step(1 - a2)* step(exp(-a0) -1 - a2)*step(1 + a1 + a2 + a3)*step(exp(- a0) - 1 - a1 - a2 - a3)} I get an error. Is this what you are suggesting? – Morris

The error message says: ...command #Bugs:set cannot be executed (is greyed out) ...Any help would be much appreciated. Thanks – Morris

Yeah, that is what I was suggesting. – Outstation

Any thoughts on how to solve this error? – Morris

It would seem that to repeat observations according to your weights you could use newData <- c(N = sum(weights), lapply(data.weighted[2:5], rep, times = weights)). – Riband

Thanks, @JuliusVainora. But the problem with that is that weights are usually real numbers, not integers. – Morris

I see, although that contradicts your "The weights variable in this [my] case aims to "[indicate] the number of times that row was observed". – Riband

That statement The weights variable in this [my] case aims to "[indicate] the number of times that row was observed is just an approximation of the truth. Likewise, newData <- c(N = sum(weights), lapply(data.weighted[2:5], rep, times = weights)) would only approximate the truth. Too much rounding errors and other associated bias. – Morris