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