multiple likelihoods with (partial) common parameters

[Scott Foster]

Dear INLA experts,

I'm wondering if this model can be estimated in INLA.  It seems simple, but may be just complex enough in just the wrong sort of way.

2 outcomes with different likelihoods.

The same set of covariates are measured at all observations.

The covariate values are different for all observations.

The linear predictor for the first likelihood is g(E(y1)=alpha1+X'beta, where alpha 1 is an intercept and beta is a vector of coefficients.

The linear predictor for the second likelihood is g(E(y2))=alpha2+gamma(X'beta), where alpha2 is another intercept, beta is as before, and gamma is a scaling coefficient.

I have made a little dummy data set (below) that highlights the problem.  Note that the real application will also require use of a random effect, shared between data types -- a spatial effect with the unique measurement locations.

Thanks in advance,

Scott

#based on the example in https://becarioprecario.bitbucket.io/inla-gitbook/ch-INLAfeatures.html

#model
#poisson model is exp( E( y1)) = alpha1 + b1*x1 + b2*x2 + b3*x3  #log link
#gaussian model is E( y2) = alpha2 + g * ( b1*x1 + b2*x2 + b3*x3) #identity link

library( INLA)

set.seed( 747)

#data attributes and covariates
N <- c(4000, 2000)
p <- 3
X <- matrix( runif( sum(N)*p,-1,1), nrow=sum(N), ncol=p, dimnames=list(NULL,c("x1","x2","x3")))
X <- model.matrix( ~-1+x1+x2+x3, as.data.frame( X))
X <- cbind( I1=1, I2=1, X) #adding another intercept
X[N[1]+1:N[2],"I1"] <- NA
X[1:N[1],"I2"] <- NA

#coefficients for data generation
coefs <- matrix( NA, ncol=2, nrow=4, dimnames= list(c("alpha","b1","b2","b3"), NULL))
coefs[,1] <- c(0.5, -1.25 ,0.25, 1)
coefs[,2] <- c( -0.25, 2.25 * coefs[-1,1])

#model products and simulated outcomes
lps <- evs <- y <- matrix( NA, nrow=sum(N), ncol=2)
#poisson data
lps[,1] <- X[,c("I1","x1","x2","x3")] %*% coefs[,1]
evs[,1] <- exp( lps[,1])
y[1:N[1],1] <- rpois( N[1], evs[1:N[1],1])
#gaussian data
lps[,2] <- X[,c("I2","x1","x2","x3")] %*% coefs[,2]
evs[,2] <- lps[,2]
y[N[1]+1:N[2],2] <- rnorm( n=N[2], sd=1, mean=evs[N[1]+1:N[2],2])

fm <- list()
#simple model with shared covariates
fm[[1]] <- inla(y ~ -1 + (I1 + I2) + (x1 + x2 + x3),
                 data = list(Y = y, I1=X[,1], I2=X[,2], x1=X[,3], x2=X[,4], x3=X[,5]),
                 family = c("poisson", "gaussian"))
summary( fm[[1]])

#simple model with distribution specific covariates (same as 2 glms - no shared parameters)
fm[[2]] <- inla(y ~ -1 + (I1 + I2) + (I1 + I2):(x1 + x2 + x3),
                data = list(Y = y, I1=X[,1], I2=X[,2], x1=X[,3], x2=X[,4], x3=X[,5]),
                family = c("poisson", "gaussian"))
summary( fm[[2]])
tmpCoef <- fm[[2]]$summary.fixed$mean[3:8]
print( tmpCoef[4:6] / tmpCoef[1:3]) #naive estimate of g

#Note that the above model does not contain a random effect.
#eventually, this will be needed (just one).

[Håvard Rue]

its doable. you need to store X'beta in another variable, v, say, using fake
data and building another linear predictor. then you can use 'v' as a part to
explain y1 and then copy v to get beta*v

I cannot recall if there is an example of this in the book...

> --
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/27e8f963-b535-4f8f-94b1-f18ed07a545fn%40googlegroups.com
> .

--
Håvard Rue
hr...@r-inla.org

[Helpdesk (Haavard Rue)]

this is an example

On Sun, 2025-05-04 at 21:11 -0700, Scott Foster wrote:

> --
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/27e8f963-b535-4f8f-94b1-f18ed07a545fn%40googlegroups.com
> .

--
Håvard Rue

he...@r-inla.org

[Scott Foster]

Dear Haavard,

Thanks for the response (and helpful code).  Apologies for my tardy reply -- hopefully there will be enough info here to jog your memory.

The code nearly does what I expect it to.  However, there are some situations where it gives unintuitive estimates.  The attached minor re-working of your example code is an illustrative example.  The problem is exacerbated when all the regression coefficients are positive.  If beta is negative (as in the example), then estimation proceeds as expected.

Some questions about the code:

* Why weight the fake zero observations with -1?  Why not '0' so that they don't contribute to the logl.  With -1, won't they (partially) cancel the contributions from the non-fake data?

* In the design matrix, the intercept term was filled with the simulation parameter.  Should this be 1 instead?  Note that Intercept.2 is assigned zero.

* The plots show that there might be a problem that the beta parameter is being applied to both otucomes.  I can't see how the model specification allows that though.

I also want to check: the covariates (x and xx in example) are the same for y1 and y2. This is not the case for my real application, nor are there the same amount of observations.  I presume that this is not a problem?

Thanks again,

Scott

[Helpdesk (Haavard Rue)]

so what you need is this one???

On Tue, 2025-06-03 at 10:36 +1000, Scott Foster wrote:
> #poisson model is exp( E( y1)) = alpha1 + b1*x11 + b2*x21 + b3*x31  #log link
> #gaussian model is E( y2) = alpha2 + g * ( b1*x12 + b2*x22 + b3*x32) #identity



--
Håvard Rue
he...@r-inla.org

[Helpdesk (Haavard Rue)]

I think this is what you want, and use different covariates as you request

On Tue, 2025-06-03 at 10:36 +1000, Scott Foster wrote:

Håvard Rue
he...@r-inla.org

[Scott Foster]

Thanks!

This is indeed what I wanted. It seems to work well.

I'm not sure why the two codes produce different results though. The first code, which is not-quite-right, copies only the data for one of the two likelihoods (the one that has to be altered). The second code, which works, copies all the data from both likelihood models. It seems to my eyes, that the model should be the same, yet different results? What am I missing?

Apart from understanding, another reason for me asking is that my 'real' data is quite numerous (possibly 10,000s of observations in a geostats model) and I wonder if it is necessary to repeat all the data, or just the relatively small section that requires this adjustment. However, I am happy enough with the computational effort if the model is the correct one.

Thanks again,

Scott

[Helpdesk (Haavard Rue)]

On Fri, 2025-06-13 at 10:38 +1000, Scott Foster wrote:
> Thanks!
>
> This is indeed what I wanted. It seems to work well.
>
> I'm not sure why the two codes produce different results though. The first
> code, which is not-quite-right, copies only the data for one of the two
> likelihoods (the one that has to be altered). The second code, which works,
> copies all the data from both likelihood models. It seems to my eyes, that the
> model should be the same, yet different results? What am I missing?

I guess this boils down to writing out the details in a simple low-dimensional
case...

>
> Apart from understanding, another reason for me asking is that my 'real' data
> is quite numerous (possibly 10,000s of observations in a geostats model) and I
> wonder if it is necessary to repeat all the data, or just the relatively small
> section that requires this adjustment. However, I am happy enough with the
> computational effort if the model is the correct one.

10,000's are ok. the 'copy' operator can only copy the whole model, but the
added cost is very small, since if w is a copy of x, then it just say that

w[i] - x[i] = small noise

so it does not add much complexity in the graph

best
H

--
Håvard Rue
he...@r-inla.org