[ADMB Users] Random effect as finite mixture
Shige Song
shigesong at gmail.com
Mon Feb 23 05:39:26 PST 2009
Dear Dave,
You are right. This can be used as a template file for handing growth
mixture modeling using ADMB. Thanks for the effort.
I plan to spend some time adapting the code to handle my empirical questions
at hand, and comparing the results to those obtained from Mplus, as well as
some benchmarks. I will report to the community once I have something to
report.
Many thanks.
Best,
Shige
On Sun, Feb 22, 2009 at 4:32 PM, dave fournier <otter at otter-rsch.com> wrote:
>
> Shige,
>
> I think the "ADMB team" consists mostly of me. So I'll try and answer
> your question about finite mixtures and growth models. Until today I had
> very little idea what growth models were, but in writing the simulator I
> hope I am getting the idea.
>
> Actually finite mixtures are handled directly so it is not necessary to
> use the random effects package. This is actually an advantage since it
> greatly speeds up the calculations.
> To illustrate the ideas I built a simulator to create a data set
> something like you describe in
> http://sgsong.blogspot.com/2009/01/finite-mixture-modeling-with-admb.html
>
> This could be done in R but I used ADMB since it illustrates some of the
> techniques which are useful for building the analyzer. since you didn't
> identify the growth curves you used I use the von Bertalanfy curve
>
> L_t = Linf * ( 1 - exp(-k*age))
>
> where L_0=0 is assumed. Of course it is a simple matter to switch to a
> different curve There are two curves one with Linf=100 and the other
> with Linf=120. For both cases k=0.11 (time measured in years).
>
> I generated 2000 samples. First the size measured in 2 month intervals
> and then at ages 8,11,15,19
>
> The individuals were generated with a probability of 0.4 of belonging to
> group1 (linf=100) and 0.6 of belonging to group 2. After the first
> period the individual were changed at random to groups 1 or 2 with a
> transition matrix
>
> .8 prob stay in group 1
> .2 prob to move from group 2 to group 1
> .7 prob to say in group 2
> .3 prob to move from group 2 to group 1
>
> The the size data was generator using the appropriate growth curve.
> the resulting sizes were then perturbed using multiplicative lognormal
> random variables (std dev of .08) which were highly positively
> correlated (0.9)
>
> The resulting output was analyzed with a finite mixture model where
> the 12 initial measuremeasurements are assumed to come from a mixture of
> two multivariate log-normal distributions with a correlation matrix of
> the form
>
> 1 rho rho^2 .... rho^11
> rho 1
> rho^2 .......
>
> etc.
>
>
> The model was fit and the Hessian evaluated in under a minute.
>
> these are the parameter estimates for a typical run.
>
> # Number of parameters = 10 Objective function value = -84436.9
> Maximum gradient component = 0.0208847
> # p1coff:
> 0.379774 0.620226
> # p2coff:
> 0.482844 0.517158
> # Linf:
> 99.2992 119.646
> # k:
> 0.109928 0.110294
> # rho:
> 0.882750993879
> # lsigma:
> -2.54165198426
>
>
> For each individual the conditional probabilities of belonging in either
> of the two groups was calculated after observing the first sizes and
> again after observing the second group of sizes. this was used to
> estimate the transition probabilities. (Of course this could be done in
> the model which would be more efficient.)
> ~
> ~ 0.793901 0.206099
> 0.292103 0.707897
> ~
> which is almost perfect.
>
> the C++ code for the simulator is
>
> ~
> ~
> #include <admodel.h>
>
> dvector growth_function(const dvector & age,const dvector& theta)
> {
> dvector tmp=theta(1)*(1.0-exp(-theta(2)*age)); // use you own
> return tmp; // growth function here
> }
>
> int main(int argc,char * argv[])
> {
> int n=2000;
> int nm=16;
> int iseed=289;
> random_number_generator rng(iseed);
> dvector p1(1,2); // priori probability that a child is in group 1 or 2
> p1(1)=.4; p1(2)=.6;
> dvector p2(1,2);
> dmatrix M(1,2,1,2); // transition probabilities
> M(1,1)=.8; // 1 remains 1
> M(2,1)=.3; // 2 moves to 1
> M(1,2)=.2; // 1 moves to 2
> M(2,2)=.7; // 2 remains 2
> dvector choose(1,n);
> dvector choose1(1,n);
> double sigma=0.08;
> double rho=0.9;
> dmatrix size(1,n,1,nm);
> imatrix group(1,n,1,2);
> dvector age(1,nm);
> dvector eta(1,nm);
> dvector eps(1,nm);
> age.fill_seqadd(.2,.2);
> age(13)=8;
> age(14)=11;
> age(15)=15;
> age(16)=19;
> dmatrix theta(1,2,1,2); //parameters for growth curves
> theta(1,1)=100;
> theta(2,1)=120;
> theta(1,2)=.11;
> theta(2,2)=.11;
>
> choose.fill_randu(rng); // use to pick group that each child belongs to
> choose1.fill_randu(rng); // use to pick group that each child belongs to
> int i;
> for (i=1;i<=n;i++) // loop over children
> {
> if (choose(i)<p1(1))
> {
> group(i,1)=1; //child in group1
> }
> else
> {
> group(i,1)=2; // child in group2
> }
> if (group(i,1)==1) // use probability that group 1 remains in group1
> if (choose1(i)<M(1,1))
> group(i,2)=1; // remained in group1
> else
> group(i,2)=2; // moved to group2
> else
> if (choose1(i)<M(2,1)) // use probability that group 2 changes to
> group1
> group(i,2)=1; // changed to group1
> else
> group(i,2)=2; // remained in group2
> }
>
> for (i=1;i<=n;i++) // loop over children and get size data
> {
> size(i)(1,12)=growth_function(age(1,12),theta(group(i,1)));
> size(i)(13,16)=growth_function(age(13,16),theta(group(i,2)));
> }
> double rho2=rho*rho;
> for (i=1;i<=n;i++) // loop over children and add noise
> {
> eta.fill_randn(rng);
> eps(1)=eta(1);
> int j;
> for (j=2;j<=4;j++) // loop over children and add noise
> {
> eps(j)=rho*eps(j-1)+sqrt(1.0-rho2)*eta(j);
> }
> eps(5)=eta(1);
> for (j=6;j<=16;j++) // loop over children and add noise
> {
> eps(j)=rho*eps(j-1)+sqrt(1.0-rho2)*eta(j);
> }
> size(i)=elem_prod(size(i),exp(sigma*eps)); // lognormal error
> }
>
> ofstream ofs("growth.dat");
> ofs << "# nobs" << endl << n << endl ;
> ofs << "# nmeasure" << endl << nm << endl ;
> ofs << "# ages" << endl << setprecision(3) << setw(7) << age << endl ;
> ofs << "# observations" << endl << setprecision(3) << setw(7) << size
> << endl ;
> ivector icount1(1,2);
> icount1.initialize();
> ivector icount2(1,2);
> icount2.initialize();
> for (i=1;i<=n;i++) // loop over children
> {
> icount1(group(i,1))++;
> icount2(group(i,2))++;
> }
> cout << "proportions in groups " << icount1/sum(icount1) << endl;
> cout << "proportions in groups " << icount2/sum(icount2) << endl;
> }
>
>
> while the TPL code for the analyzer is
>
> DATA_SECTION
> init_int n
> init_int nm
> init_vector ages(1,nm)
> init_matrix sizes(1,n,1,nm)
> matrix q0(1,n,1,2)
> matrix q1(1,n,1,2)
> number lmin
> number lmax
> LOC_CALCS
> ages.fill_seqadd(.2,.2);
> ages(13)=8;
> ages(14)=11;
> ages(15)=15;
> ages(16)=19;
> lmin=min(column(sizes,nm));
> lmax=max(column(sizes,nm));
> PARAMETER_SECTION
> init_bounded_vector p1coff(1,2,.001,1.1)
> init_bounded_vector p2coff(1,2,.001,1.1)
> vector p1(1,2)
> vector p2(1,2)
> number psum1
> number psum2
> init_bounded_vector Linf(1,2,0.9*lmin,1.1*lmax,2)
> init_bounded_vector k(1,2,.01,.3)
> init_bounded_number rho(-.1,.95,3)
> init_bounded_number lsigma(-5.0,5.0)
> number sigma
> matrix cor1(1,12,1,12)
> matrix cor2(1,4,1,4)
> matrix is1(1,12,1,12)
> matrix is2(1,4,1,4)
> matrix ic1(1,12,1,12)
> matrix ic2(1,4,1,4)
> matrix pred_size(1,2,1,nm)
> number lds1
> number lds2
> objective_function_value f
> PROCEDURE_SECTION
> get_correlation();
> get_proportions();
> get_predicted_sizes();
> get_likelihood();
>
> FUNCTION get_correlation
> f+=square(rho);
>
> sigma=exp(lsigma);
> dvariable vinv=1.0/(sigma*sigma);
> int i,j;
> for (i=1;i<=12;i++)
> for (j=1;j<=12;j++)
> cor1(i,j)=pow(rho,fabs(i-j));
>
> for (i=1;i<=4;i++)
> for (j=1;j<=4;j++)
> cor2(i,j)=pow(rho,fabs(i-j));
>
> ic1=inv(cor1);
> ic2=inv(cor2);
>
> is1=ic1*vinv;
> is2=ic2*vinv;
>
> lds1=ln_det(is1);
> lds2=ln_det(is2);
>
> FUNCTION get_proportions
> psum1=sum(p1coff);
> p1=p1coff/psum1;
> f+=100.*square(log(psum1));
>
>
> psum2=sum(p2coff);
> p2=p2coff/psum2;
> f+=100.*square(log(psum2));
>
> FUNCTION get_predicted_sizes
> lmax=max(column(sizes,nm));
> pred_size(1)=log(Linf(1)*(1.0-exp(-k(1)*ages)));
> pred_size(2)=log(Linf(2)*(1.0-exp(-k(2)*ages)));
>
> FUNCTION get_likelihood
> int i;
> dvar_matrix res(1,2,1,nm);
> for (i=1;i<=n;i++)
> {
> res(1)=log(sizes(i))-pred_size(1);
> res(2)=log(sizes(i))-pred_size(2);
> dvar_vector res1a=res(1)(1,12);
> dvar_vector res1b=(res(1)(13,16)).shift(1);
> dvar_vector res2a=res(2)(1,12);
> dvar_vector res2b=(res(2)(13,16)).shift(1);
> dvariable l1a=exp(-0.5*res1a*(is1*res1a));
> dvariable l1b=exp(-0.5*res1b*(is2*res1b));
> dvariable l2a=exp(-0.5*res2a*(is1*res2a));
> dvariable l2b=exp(-0.5*res2b*(is2*res2b));
> f-=log(1.e-10+(p1(1)*l1a+p1(2)*l2a));
> f-=log(1.e-10+(p2(1)*l1b+p2(2)*l2b));
> f-=0.5*(lds1+lds2);
>
> q0(i,1)=value(p1(1))*value(l1a);
> q0(i,2)=value(p1(2))*value(l2a);
> q0(i)/=sum(q0(i));
> q1(i,1)=value(p2(1))*value(l1b);
> q1(i,2)=value(p2(2))*value(l2b);
> q1(i)/=sum(q1(i));
> }
> REPORT_SECTION
> report << setprecision(4) << setw(8) << inv(is1) << endl;
> report << setprecision(4) << setw(8) << inv(is2) << endl;
>
> report << setprecision(4) << setw(8) << cor1 << endl;
> report << setprecision(4) << setw(8) << cor2 << endl;
>
> report << colsum(q0)/n << endl;
> report << colsum(q1)/n << endl;
>
> report << n << endl;
> report << q0 << endl;
> report << q1 << endl;
>
>
> I hope this illustrates how ADMB might be applied to growth models.
>
> Dave
>
>
>
>
>
> --
> David A. Fournier
> P.O. Box 2040,
> Sidney, B.C. V8l 3S3
> Canada
> Phone/FAX 250-655-3364
> http://otter-rsch.com
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20090223/49a00b13/attachment.html>
More information about the Users
mailing list