Richard Chandler richard.chandlers at gmail.com
Fri Mar 11 07:40:08 PST 2011

```Thank you Hans and Dave. The tpl file below now works. Weihai Liu was
also kind enough to send me a suggestion, which is pasted in too.
Since Hans is one of the principal developers, I assume his
implementation is the preferred one.

Richard

// Based on Hans' suggestions

PARAMETER_SECTION

init_bounded_number log_lambda(-20.0,20.0,1)
init_bounded_number p0(-20.0,20.0,1)
init_bounded_number p1(-20.0,20.0,1)

init_bounded_number log_sigma(-3.0,3.0,2) // log of random effect SD
random_effects_vector u(1,nG,2)

objective_function_value nll

PROCEDURE_SECTION

for(int id=1; id<=nG; id++) {
nll_group(id, p0,p1,log_lambda,log_sigma,u(id));
}

SEPARABLE_FUNCTION void nll_group(int id, const dvariable& p0,const
dvariable& p1,const dvariable& log_lambda, const dvariable& log_sigma,
const dvariable& u)
dvar_vector p(1,R);
dvar_vector f(1,S);
dvar_vector g(1,S);
dvariable sigma = exp(log_sigma);
dvariable lambda = exp(log_lambda);

nll += log_sigma + 0.5*square(u/sigma);

int q=10*(id-1)+1;
for(int i=q;i<=q+9;i++) {
p(i) = 1.0/(1.0+exp(-1.0*(p0 + p1*x(i) + u)));
for(int k=1;k<=S;k++) {
f(k) = pow(lambda, N(k)) * exp(-1.0*lambda - gammln(N(k)+1));
g(k) = 1.0;
for(int j=1;j<=T;j++) {
if(N(k)>=y(i,j))
g(k) *= exp(gammln(N(k)+1) - gammln(y(i,j)+1) -
gammln(N(k)-y(i,j)+1)) * pow(p(i), y(i,j)) *
pow(1.0-p(i), N(k)-y(i,j));
else
g(k) =0.0;
}
}
nll -= log(f*g);
}

// based on Weihai's example

PARAMETER_SECTION

init_bounded_number log_lambda(-20.0,20.0,1)
init_bounded_number p0(-20.0,20.0,1)
init_bounded_number p1(-20.0,20.0,1)

init_bounded_number log_sigma(-3.0,3.0,2) // log of random effect SD
random_effects_vector u(1,nG,2)

objective_function_value nll

PROCEDURE_SECTION
nll_group(p0,p1,log_lambda,log_sigma,u);

SEPARABLE_FUNCTION void nll_group(const dvariable& p0,const dvariable&
p1,const dvariable& log_lambda, const dvariable& log_sigma, const
dvar_vector& u)
dvar_vector p(1,R);
dvar_vector f(1,S);
dvar_vector g(1,S);
dvariable sigma = exp(log_sigma);
dvariable lambda = exp(log_lambda);

//  nll = nG*log_sigma + 0.5*norm2(u/sigma);
nll = nG*log_sigma + 0.5*norm2(u)/mfexp(2*log_sigma);

for(int i=1;i<=R;i++) {
int id = ID(i);
p(i) = 1.0/(1.0+exp(-1.0*(p0 + p1*x(i) + u(id))));
for(int k=1;k<=S;k++) {
f(k) = pow(lambda, N(k)) * exp(-1.0*lambda -
gammln(N(k)+1)); //poisson
g(k) = 1.0;
for(int j=1;j<=T;j++) {
if(N(k)>=y(i,j))
g(k) *= exp(gammln(N(k)+1) - gammln(y(i,j)+1) -
gammln(N(k)-y(i,j)+1)) * pow(p(i), y(i,j)) *
pow(1.0-p(i), N(k)-y(i,j));//binomial
else
g(k) =0.;
}
}
nll -= log(f*g + 1.e-20);
}

On Fri, Mar 11, 2011 at 6:08 AM, H. Skaug <hskaug at gmail.com> wrote:
> Hi,
>
> Yes, it seems that your model is separable, but of type "4.1 The first example"
> in ADMBRE manual. Hence, you do not need the sparse matrix stuff.
> In your code you need to loop over the elements of "u" like in "4.1
> The first example"
> Your PROCEDURE_SECTION should basically consist of:
>
>  for(int id=1;id<=xxxxx;id++)
>   SEPARABLE_FUNCTION(u(id), + the rest of your arguments)
>
>
> Inside the SEPARABLE_FUNCTION you need to do all the calculations
> that involve u(id).
>
> If you do this you will be able to handle large models. Make sure
> that you get the same result on small examples with and with SEPARABLE_FUNCTION.
>
>
> Hans
>
> On Thu, Mar 10, 2011 at 2:54 PM, Richard Chandler
> <richard.chandlers at gmail.com> wrote:
>> Hi,
>>
>> I implemented a binomial-Poisson mixture model (Royle Biometrics 2004)
>> with random effects in ADMB (.tpl below). It works very well on small
>> simulated datasets; however, I run into memory issues when I try
>> complex variations of the model on large datasets. I have fiddled with
>> many of the command line options and the TOP_OF_MAIN_SECTION, but I
>> believe I need to use ADMB's sparse matrix capabilities via the
>> simple:
>>
>> SEPARABLE_FUNCTION void nll_group(const dvariable& log_sigma, const
>> dvar_vector& u)
>>  nll = nG*log_sigma + 0.5*norm2(u)/mfexp(2*log_sigma);
>>
>> but this does not converge. I am new to ADMB and would greatly
>> appreciate it if someone could suggest a better way to use a
>> SEPARABLE_FUNCTION here. The RE manual is very helpful, but I must be
>> missing something. Also, I would appreciate any other suggestions for
>> improving the implementation of this model.
>>
>> Thanks for the great software.
>>
>> Richard
>>
>

```