dave fournier davef at otter-rsch.com
Fri May 9 13:27:37 PDT 2014

```I noticed that someone had posted  on the SAS list about a moderate
sized model
requiring something like 80 hours to fit using GLIMMIX.  I thought this
might be
an interesting application of some of the sparse matrix stuff.  The
posted a simple example of the kind of problem he was having trouble with.

http://listserv.uga.edu/cgi-bin/wa?A2=ind1405b&L=sas-l&P=4001

He also later said that his real problem has an ar(1) structure which we
know how to deal with.
(I think).  I wrote the attached TPL file to test the ideas.

I ran the model with the following.

time ./ryan -mno 1000000 -ams 10000000 -ilmn 5 -shess -ndi 100000
-noinit -l1 100000000 -nl1 100000000  -iprint 1  -crit 1.e-3 -noreders

the -noreders  is to prevent the calculation of the estimates of the std
devs of the random effects via the delta method.

This saves a lot of time and space since there about 75,000 random
effects.  Fitting the model required less than 4 minutes.
I don't know how many ADMB people are interested in this kind of
larger model.

It appears that the -noreders option has never been included in the
code.  It is only about 8 or so lines of code to add.
Maybe someone is interested.  Per usual I will hold my breath.
-------------- next part --------------
DATA_SECTION
int nloc
int nsub
int nobs
int iseed
LOC_CALCS
nloc=250;
nsub=300;
nobs=4;
END_CALCS
3darray y(1,nloc,1,nsub,1,nobs)
matrix x(1,nloc,1,nsub)
LOC_CALCS
iseed=12345;
random_number_generator rng(iseed);

double lre;
double sre;
double prob;
double eta;

for (int i=1;i<=nloc;i++)
{
lre=sqrt(.3)*randn(rng);
for (int j=1;j<=nsub;j++)
{
sre=sqrt(.5)*randn(rng);
x(i,j)=randn(rng);
eta=1.6+2.5*x(i,j)+lre+sre;
prob=exp(eta)/(1.0+exp(eta));
for (int k=1;k<=nobs;k++)
{
double z=randu(rng);
if (z<=prob)
y(i,j,k)=1.0;
else
y(i,j,k)=0.0;
}
}
}

PARAMETER_SECTION
objective_function_value f
init_number a
init_number b
init_bounded_number lsd(0.0,5.0,2)
init_bounded_number ssd(0.0,5.0,2)
!! lsd.set_scalefactor(100.);
!! ssd.set_scalefactor(100.);
random_effects_vector eloc(1,nloc,2)
random_effects_matrix esub(1,nloc,1,nsub,2)

PROCEDURE_SECTION

for (int i=1;i<=nloc;i++)
{
for (int j=1;j<=nsub;j++)
{
f1(i,j,a,b,lsd,ssd,eloc(i),esub(i,j));
}
}

SEPARABLE_FUNCTION void f1(int i,int j,const prevariable& a,const prevariable& b,const prevariable& lsd,const prevariable& ssd,const prevariable& el,const prevariable& es)

dvariable eta=a+b*x(i,j)+lsd*el+ssd*es;
dvariable prob;
if (value(prob)<20.)
prob=exp(eta)/(1.0+exp(eta));
else
prob=1.0/(1.0+exp(-eta));

for (int k=1;k<=nobs;k++)
{
if (y(i,j,k)>0)
f-=log(1.e-10+prob);
else
f-=log(1.e-10+(1.0-prob));
}
if (j==1)
{
f+=0.5*square(el);
}
f+=0.5*square(es);

```