Interactions and linear combinations in a joint likelihood model

[Tyler Scott]

Greetings all,

    

I have a question about the proper procedure for examining interactions within the context of a joint likelihood model. The following is example code for a model where the first stage is a 0/1 binary dependent variable, and then the second stage with a continuous dependent variable only for observations where u_i == 1 at the first stage. The basic procedure borrows from the code found in the Blangiardi/Cameletti book. In any case, I wish to generate and then plot marginal effects for an interaction term that is present in both stages of the model. The way that I can get the models with linear combinations (rr and rr2) to run without an error message is to index the linear combinations row-wise by stage (i.e., 1:n, or (n+1:n*2)). 

So, I have two (probably rudimentary) questions:

1. Is this procedure correct, or is something different required for the joint likelihood model? 

n = 100

u <- (rnorm(n)>0) + 0

y <- rgamma(n,10)

y[u==0] <- NA

idat <- list(Y=matrix(NA,2*n,2))

idat$Y[1:n,1] <- u

idat$Y[n+1:n,2] <- y

idat$mu.u <- rep(1:0, each=n)

idat$mu.y <- rep(0:1, each=n)

idat$u_violation <- c((rnorm(n)>0) + 0, rep(0,n))

idat$y_violation <- c(rep(0,n),idat$u_violation[1:n])

idat$u_revenue <- c(runif(n,0,100), rep(0,n))

idat$y_revenue <- c(rep(0,n),idat$u_revenue[1:n])

form = Y ~ 0 + mu.u + mu.y + u_revenue + y_revenue + u_violation + y_violation + 

  u_revenue:u_violation + y_revenue:y_violation

r = inla(form,c('binomial', 'gamma'),

           data=idat,control.predictor=list(compute=TRUE),verbose=F)

lc_u = inla.make.lincombs( 

  u_revenue =  r$model.matrix[1:n, "u_revenue"], 

  `u_revenue:u_violation` =  r$model.matrix[1:n, "u_revenue:u_violation"]) 

lc_y = inla.make.lincombs( 

  y_revenue =  r$model.matrix[(n+1):(n*2), "y_revenue"], 

  `y_revenue:y_violation` =  r$model.matrix[(n+1):(n*2), "y_revenue:y_violation"]) 

rr = inla(form,c('binomial', 'gamma'),

    data=idat,control.predictor=list(compute=TRUE),verbose=F,lincomb = lc_u)

rr2 = inla(form,c('binomial', 'gamma'),

          data=idat,control.predictor=list(compute=TRUE),verbose=F,lincomb = lc_y)

2. Am I correct in assuming then that the following code then plots the posterior density of the revenue parameter for the first stage of the model, for non-violators and then for violators?

plot(density(rr$summary.lincomb.derived$mean[idat$u_violation[1:100]==0]),ylim=c(0,2))

lines(density(rr$summary.lincomb.derived$mean[idat$u_violation[1:100]==1]),col='red')

thank you for your time.

[INLA help]

Hi

this seems correct. You can, if you want, merge, lc_y and lc_u, if you
give them different names, so you only need one call to inla() again,
like

> names(lc_u) = paste('u', 1:length(lc_u), sep="")
> names(lc_y) = paste('y', 1:length(lc_y), sep="")
> lc = c(lc_u, lc_y)


but if you're unsure about the valididity of the results, just compute
the mean using the mean of fixed effects, and, artificially just make
sure the interaction is orthogonal or almost so, to the fixed effects
(just for testing) and then the variance is just the sum.

> 2. Am I correct in assuming then that the following code then plots
> the posterior density of the revenue parameter for the first stage of
> the model, for non-violators and then for violators?
>
> plot(density(rr$summary.lincomb.derived$mean[idat$u_violation[1:100]=
> =0]),ylim=c(0,2))
> lines(density(rr$summary.lincomb.derived$mean[idat$u_violation[1:100]
> ==1]),col='red')

this I am unsure about, since

> sum(idat$u_violation[1:100]==0)
[1] 53


but the plot, will only plot one marginal. may you want a loop ???


best
H


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

[Tyler Scott]

Hello,

    Thanks for your prior response. Back with a related issue. Using the same sample code:

n = 100; u <- (rnorm(n)>0) + 0; y <- rgamma(n,10); y[u==0] <- NA; dat <- list(Y=matrix(NA,2*n,2))

dat$Y[1:n,1] <- u

dat$Y[n+1:n,2] <- y

dat$mu.u <- rep(1:0, each=n)

dat$mu.y <- rep(0:1, each=n)

violation = (rnorm(n)>0) + 0

dat$u_violation <- c(violation,rep(NA,n)); dat$y_violation <- c(rep(NA,n),violation)

dat$u_revenue <- c(runif(n,0,100), rep(0,n)); dat$y_revenue <- c(rep(0,n),dat$u_revenue[1:n])

r = inla(Y ~ 0 + mu.u + mu.y + u_revenue*u_violation + y_revenue*y_violation, c('binomial', 'gamma'),data = dat, control.predictor = list(compute=TRUE)) 

lc = inla.make.lincombs(u_violation= r$model.matrix[1:n, "u_violation"],  "u_revenue:u_violation" = r$model.matrix[1:n, "u_revenue:u_violation"]) 

rr = inla(Y ~  0 + mu.u + mu.y +  u_revenue*u_violation + y_revenue*y_violation,c('binomial', 'gamma'), data = dat , control.predictor = list(compute=TRUE), lincomb = lc) 

I now wish to compute a posterior density for `violation + violation:revenue`, where violation is a dummy variable (basic idea is the effect of a violation on the response variable varies by revenue level). In any case, the code above produces a crash at the first linear combination where both violation and violation*revenue are 0 (i.e., the first non-violator, which makes sense since the marginal effect is then zero as well).

*** ERROR ***	Lincomb error for section[lincomb.lc001]
*** ERROR ***	All weights are zero.
	GMRFLib version 3.0-0-snapshot, has recived error no [9]
	Reason    : Invalid value of parameter
	Message   : Condition `0 == 1' is not TRUE
	Function  : inla_parse_lincomb
	File      : src/inla.c
	Line      : 8948
	RCSId     : hgid: 7e45fc2b7a2f  date: Tue Jan 31 09:14:12 2017 +0300
Aborted
Error in inla.inlaprogram.has.crashed() : 
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <he...@r-inla.org>.

My question then is this: What is the appropriate procedure in this case? Is it appropriate to simply compute linear combinations only for violators (as in the code below), or does that raise other problems?

lc_violators = inla.make.lincombs(
  u_violation= r$model.matrix[which(violation==1), "u_violation"], 
  "u_revenue:u_violation" = r$model.matrix[which(violation==1), "u_revenue:u_violation"]) 

rr = inla(Y ~  0 + mu.u + mu.y +  u_revenue*u_violation + y_revenue*y_violation,c('binomial', 'gamma'),
          data = dat ,
          control.predictor = list(compute=TRUE), 
          lincomb = lc_violators)

Best,

-tyler

[INLA help]

in this case you know the answer, its zero with prob=1. ;-)

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

[Tyler Scott]

Silly me, thanks for your patience Havard!