Implementing weighted log-likelihoods

[Seong Hoon Yoon]

Hi,

I am using TMB to implement optimization for my model which has random effects.

I need to attach the sampling weights to the individual log-likelihood after the log-likelihood has integrated out the random effects.

For example, I am currently using something like below to attach the weights

for(int i=0; i<N; i++){

jnll -=    weight(i) * loglik(i)

But I am wondering whether that is attaching the weights after integral approximation or before integral approximation. Random effects are specified in MakeADFun function in R.  If latter is the case,  how should I change my code so that it is attached after solving the integrals?

I highly appreciate any help.

Thank you. 

[Hans Skaug]

Your code, if put in the .cpp file, will apply the weights "before" the Laplace approximation. 

I do not think it is easy to apply the weights "after" the Laplace approximation in TMB.

On ad-hoc, and computationally inefficient, way is to call C++ once for each unique data value,

and the collect weight(i) * loglik(i) in R. To avoid overhead related to calling MakeADFun multiple times

I imagine you need to use DATA_UPDATE. This approach may not be applicable to any situation.

Hans

[Seong Hoon Yoon]

Dear Hans,

Thanks for your reply.  I have some additional questions regarding the using TMB for my model. 

1) I  have two random effects and thus in my model, there are two integrals to approximate for 'each' observation.   

However, when I specify the joint likelihood and the distributions of random effect vectors in C++ it looks like the integral approximation is done for the joint likelihood rather than individual likelihood.

So in order to implement integral approximation per observation,  I need to call C++ per observation. Is this correct?

2) Is it possible in R to collect log-likelihoods separately by calling C++ per observation (Using multiple MakADFun function or using DATA_UPDATE),  sum them, and optimize it using optim function?

Best regards,

Seong Hoon

[Hans Skaug]

In order to clarify my idea I have written an example (see attahced).

nll_weighted() is now a weighted nll that can be fed to an R function minimizer.

I have not done any serious testing of the code, so it should be taken as an "idea".

Hans

[S Yoon]

Thank you very much,  that helped a lot. I now have a better understanding of how TMB works.

[S Yoon]

Dear Hans,

I have two more questions.  I have tried your method and it works it seems to  take a lot of time.  

1. I am wondering whether it makes sense to  implement Laplace approximation in C++ manually  and  attach the weights inside C++. I have two random effect vectors, so I will need to approximate two integrals per subject. (e.g  apply Laplace approximation twice per subject)  So for N=500,  I will need to apply Laplace approximation for 1000 integrals. 

 Could this be done using a manual Laplace approximation with a for loop inside C++,   using something like https://kaskr.github.io/adcomp/laplace_8cpp-example.html ?

2. Most importantly,  I need to extract individual gradient functions and Hessians after Laplace approximation, for variance estimation.  Will there be a way to obtain these in C++ if the method above is successful?

Best regards,

Seong Hoon

[Kasper Kristensen]

Running Laplace on the C++ side has become significantly simpler with recent versions of TMB. The example https://kaskr.github.io/adcomp/TMBad_2spa_gauss_8cpp-example.html is easily modified to do the same as Hans' example - see attached.

[S Yoon]

Dear Kasper,

Thank you very much for your help. I see that Laplace implementation in C++ is simpler than expected. I am still wondering whether there is a way to extract individual gradient and Hessian functions (eg first and second derivatives of individual log-likelihoods) after manual Laplace approximation in C++, using the attached example.

I very much appreciate your help. 

[Kasper Kristensen]

As you know (from the example you linked to) there's a namespace 'autodiff' that can be used on C++ side?

In short, you could modify my example by adding a struct like this:

struct model_evaluator {
  model obj;
  ad operator()(vector<ad> sigma) {
    obj.sigma = sigma;
    return obj.eval_nldens<ad>();
  }
};

And use the struct inside the objective function like this:

    model_evaluator F = {obj};
    vector<Type> g = autodiff::gradient(F, sigma);
    matrix<Type> h = autodiff::hessian(F, sigma);

Of course, this adds overhead, so one should probably do it conditionally on some flag... And perhaps use the fact that the AD grad and hessian tapes are almost identical (no need to repeat 1000 times)... - I'll leave that to you.

[S Yoon]

Dear Kasper,

Thank you very much for your assistance.  May I ask one last question?

I am dealing with two random effect vectors which means Laplace approximation should be done two times for each individual.

In the example C++ code you attached,  I am wondering how the below part of the code should be changed to allow two random effect vectors.

For example, if I want to add another random effect vector called "u2", would it be sufficient to just add another vector and state its distribution?  My change to the code is shown in bold.

// See 'spa_gauss.cpp' 

typedef TMBad::ad_aug ad; 

struct model { 

 ad x; // Data

 vector<ad> sigma; // Parameters 

 // joint nll as function of random effects 

 ad operator()(vector<ad> u) (vector<ad> u2) { 

 ad nll = -dnorm(u,ad(0),ad(1),true).sum(); 

 ad nll = -dnorm(u2,ad(0),ad(1),true).sum(); 

 nll -= dnorm(x,sum(u),sigma(0),true); 

 return nll; 

[Kasper Kristensen]

Unfortunately that wouldn't work. The Laplace interface expects a functor with a single vector input. So, you would have to pass one vector, then split it manually within the function.

[S Yoon]

Thank you very much. I appreciate it. I think I now have much better understanding how it works.

To clarify my understanding, so If I pass a single vector of say, length 100, where the 1st to 50th element is the first random effect vector and 51st to 100th is the second one, 

would something like following work?

// See 'spa_gauss.cpp' 

typedef TMBad::ad_aug ad; 

struct model { 

 ad x; // Data

 vector<ad> sigma; // Parameters 

 // joint nll as function of random effects 

 ad operator()(vector<ad> u) { 

vector<ad> u1 = u(0:49)

vector<ad> u2 = u(50:99)

 ad nll = -dnorm(u1,ad(0),ad(1),true).sum(); 

 ad nll = -dnorm(u2,ad(0),ad(1),true).sum(); 

 nll -= dnorm(x,sum(u),sigma(0),true); 

 return nll;