[ADMB Users] promoting the admb software

[Anders Nielsen]

Hi Dave, 

I agree that the random effects part of ADMB is quite unique. At the 
course in Santa Barbara a week from now we will make sure to show how 
great the random effects part is. It is a two days introductory course,
so I think that example is a little too advanced to show there, but we 
should keep it in mind for later courses. 

I have not followed the discussion on the R-list, but am convinced that 
R has nothing with similar generality for random effect. You could post 
the example on the R-list, but when we have tried that in the past we 
have not always gotten a fair response (I remember some BS response 
about an outdated spam filter and what that shows about being careful). 
I think this example is important enough to publish in for instance 
'journal of statistical software'. I think there you would get a 
chance to reach approximately the same correct audience, it would be 
peer review so it would carry more weight.  What do you think about 
that?

Cheers,

Anders. 

 


On Sat, 2009-02-28 at 16:29 -0800, dave fournier wrote:
> Oops I sent that unfinished by mistake.
> It seems to me that r has only rudimentry support for nonlinear mixed
> models and that there are many r users with needs that could be filled
> with ADMB. See for example the recent discussion on R help aboutthe grat
> difficulty (apparently)  in incorporating random effects into
> the type of beta regression model as described in Smithson and Verkilon.
> I built a proof of concept version with nested data and
> random effects on both parameters (mean and dispersion), mu and phi in
> the language of the above paper. It sems to work well without even
> bothering to get good initial estimates for the parameters.
> I built a simulator to verify this.
> The code for the analyzer is
> 
> 
> DATA_SECTION
>   init_int nobs
>   init_int m
>   init_matrix obs(1,nobs,0,m+1)
>   int nplots
>  !! nplots=max(column(obs,m+1));
>   ivector iswitch(1,nplots)
> PARAMETER_SECTION
>   objective_function_value f
>   init_vector bet(1,m)
>   init_vector delta(1,m)
>   init_bounded_number sigma1(.01,2.0,2)
>   init_bounded_number sigma2(.01,2.0,2)
>   random_effects_vector u(1,2*nplots,2)
>   //random_effects_vector v(1,nplots,2)
> PROCEDURE_SECTION
>   int i;
>   iswitch.initialize();
>   for (i=1;i<=nobs;i++)
>   {
>     int offset=2*(obs(i,m+1)-1)+1;
>     f1(i,bet,delta,u(offset),u(offset+1),sigma1,sigma2);
>   }
> 
> SEPARABLE_FUNCTION void f1(int i,const dvar_vector& bet,const
> dvar_vector& delta, const prevariable& ui,const prevariable& vi,const
> prevariable& sigma1,const prevariable& sigma2)
> 
>     dvariable mu;
>     dvariable phi;
>     double y=obs(i,0);
>     dvariable bx=bet*obs(i)(1,m)+sigma1*ui;
>     if (value(bx)<10)
>     {
>       dvariable ebx=exp(bx);
>       mu=ebx/(1.0+ebx);
>     }
>     else
>     {
>       dvariable ebx=exp(-bx);
>       mu=1.0/(1.0+ebx);
>     }
>     phi=exp(-delta*obs(i)(1,m)+sigma2*vi);
>     f-=ln_beta_density(y,mu,phi);
>     if (iswitch(obs(i,m+1))==0)  // so that the prior only gets added once
>     {
>       f+=0.5*square(ui);
>       f+=0.5*square(vi);
>       iswitch(obs(i,m+1))=1;
>     }
> 
> REPORT_SECTION
>   //report << sqrt(norm2(u)/u.indexmax())*sigma1 << endl;
>   //report << sqrt(norm2(v)/v.indexmax())*sigma2 << endl;
> 
> GLOBALS_SECTION
> 
>   #include <admodel.h>
>   #include <df1b2fun.h>
> 
>   dvariable betaln(const prevariable& a,const prevariable& b )
>   {
>     return gammln(a)+gammln(b)-gammln(a+b);
>   }
> 
> 
>   dvariable ln_beta_density(double y,const prevariable & mu,
>     const prevariable& phi)
>   {
>     dvariable omega=mu*phi;
>     dvariable tau=phi-mu*phi;
>     dvariable lb=betaln(omega,tau);
>     dvariable d=(omega-1)*log(y)+(tau-1)*log(1.0-y)-lb;
>     return d;
>   }
> 
>   df1b2variable betaln(const df1b2variable& a,const df1b2variable& b )
>   {
>     return gammln(a)+gammln(b)-gammln(a+b);
>   }
> 
>   df1b2variable ln_beta_density(double y,const df1b2variable & mu,
>     const df1b2variable& phi)
>   {
>     df1b2variable omega=mu*phi;
>     df1b2variable tau=phi-mu*phi;
>     df1b2variable lb=betaln(omega,tau);
>     df1b2variable d=(omega-1)*log(y)+(tau-1)*log(1.0-y)-lb;
>     return d;
>   }
> 
> and the code for the simulator is
> 
> DATA_SECTION
>   init_int nplots
>   init_int nreps
>   init_int m
>   int nobs
>   init_vector bet(1,m)
>   init_vector delta(1,m)
>   init_number sigma1
>   init_number sigma2
>  !! nobs=nplots*nreps;
>   matrix obs(1,nobs,0,m+1)
>   vector u(1,nplots)
>   vector v(1,nplots)
>   matrix w(1,nplots,1,nreps)
>  LOC_CALCS
>   random_number_generator rng(921);
>   u.fill_randn(rng);
>   v.fill_randn(rng);
>   w.fill_randu(rng);
> 
>   int i=0;
>   for (int ii=1;ii<=nplots;ii++)
>   {
>     double mu;
>     double phi;
>     for (int j=1;j<=nreps;j++)
>     {
>       i++;
>       obs(i)(2,m).fill_randu(rng); // get some radnom observations
>       obs(i,1)=1; // for the mean
>       obs(i,m+1)=ii; // for the mean
>       double bx=bet*obs(i)(1,m)+sigma1*u(ii);
>       if (bx<10)
>       {
>         double ebx=exp(bx);
>         mu=ebx/(1.0+ebx);
>       }
>       else
>       {
>         double ebx=exp(-bx);
>         mu=1.0/(1.0+ebx);
>       }
>       phi=exp(-delta*obs(i)(1,m)-sigma2*v(ii));
>       double omega=mu*phi;
>       double tau=phi-mu*phi;
>       obs(i,0)=inv_cumd_beta_stable(omega,tau,w(ii,j),1.e-8);
>     }
>   }
>   ofstream ofs("output.dat");
>   ofs << "#n    m" << endl;
>   ofs << " " << nobs << "  " << m << endl;
>   ofs << "# obs " << endl;
>   ofs << obs << endl;
> 
>   ad_exit(1);
> PARAMETER_SECTION
>   objective_function_value f
>   init_vector bet(1,m)
>   init_vector delta(1,m)
> PROCEDURE_SECTION
>   int i;
>   dvariable mu;
>   dvariable phi;
>   for (i=1;i<=nobs;i++)
>   {
>     double y=obs(i,0);
>     dvariable bx=bet*obs(i)(1,m);
>     if (value(bx)<10)
>     {
>       dvariable ebx=exp(bx);
>       mu=ebx/(1.0+ebx);
>     }
>     else
>     {
>       dvariable ebx=exp(-bx);
>       mu=1.0/(1.0+ebx);
>     }
>     phi=exp(-delta*obs(i)(1,m));
>     f-=ln_beta_density(y,mu,phi);
>   }
> 
> GLOBALS_SECTION
> 
>   #include <admodel.h>
> 
>   dvariable betaln(const prevariable& a,const prevariable& b )
>   {
>     return gammln(a)+gammln(b)-gammln(a+b);
>   }
> 
> 
>   dvariable ln_beta_density(double y,const prevariable & mu,
>     const prevariable& phi)
>   {
>     dvariable omega=mu*phi;
>     dvariable tau=phi-mu*phi;
>     dvariable lb=betaln(omega,tau);
>     dvariable d=(omega-1)*log(y)+(tau-1)*log(1.0-y)-lb;
>     return d;
>   }
> 
> Finally input file for the simulator is
> 
> # nplots  nreps  m
>    50     25    2
> # bet
>    3.0  -5.0
> #delt
>   -4.0  3.0
> # sigma1  sigma2
>   .2       .5
> 
> and The parameter estimates from the model were
> # Number of parameters = 6  Objective function value = -1303.87  Maximum
> gradient component = 6.49673e-05
> # bet:
>  2.92653 -4.99330
> # delta:
>  -3.99151 2.94916
> # sigma1:
> 0.181131822789
> # sigma2:
> 0.522257834073
> # u:
> 
> The model first does a beta regression without the random effects
> in the first phase. these are the outputs One can see that adding the
> random effects  improved the parameter estimates considerably and also
> led to a big increase in the log-likelihood
> 
> # Number of parameters = 4  Objective function value = -1136.36  Maximum
> gradient component = 2.75648e-06
> # bet:
>  2.91619 -4.96991
> # delta:
>  -3.62180 2.67383
> 
> 
> 
> 
> 
> 
> 
> 
> 
>