[ADMB Users] promoting the admb software

[dave fournier]

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










-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com