I am trying to implement a non linear optimisation for a non-hidden multi-state markov model. A similar methodology, if you are interested, can be found in the msm R package (methodology explained in https://cran.r-project.org/web/packages/msm/vignettes/msm-manual.pdf).
the input data are composed by:
TIMELAG NOCC TR1 TR2 TR3 TR4
timelag is a float value indicating the timelag of the transition
nocc is an integer indicating how many times the transition is happening
tr1 is 0/1 depending on whether this transition is 0-0 or not
tr2 is 0/1 depending on whether this transition is 0-1 or not
tr3 is 0/1 depending on whether this transition is 1-0 or not
tr4 is 0/1 depending on whether this transition is 1-1 or not
(I need the 4 transition variables to understand with IML which record to select of the Pmatrix below)
what I want to optimise is basically a piecewise likelihood function, for every record I will have a contribution, then I'll sum all contributions and that will be my likelihood function.
I start with take 2 values as input, here e.g.:
0.3, 0.7
then compose the Qmatrix as:
-0.3 0.3 ,
0.7 -0.7
the create the Pmatrix for each record as:
exp(timelag*qmat11) exp(timelag*qmat11)*qmat12
exp(timelag*qmat22)*qmat21 exp(timelag*qmat22)
then depending on which type of transition it is, I use as contribution to the likelihood the respective record in the pmatrix, e.g. if 0-0 I select Pmat11.
since this computation has to happen for every record, I thought it might have been much simpler to do that rowwise, so instead of using 2x2 matrixes, I will have a vector with 4 elements as qmatrix for each record, a vector with 4 elements as pmatrix for each record, then sum the product between pmatrix columns and tr1 tr2 tr3 tr4 columns, in order to save only the correct contribution.
proc iml;
reset noname;
use mod10;
read all var {WHICHA TIMELAG NOCC} into myds;
read all var {TR1 TR2 TR3 TR4} into trds;
close mod10;
start mylik(x) global(myds,trds);
qmat=shape(-x[1],nrow(myds),1)||shape(x[1],nrow(myds),1)||shape(x[2],nrow(myds),1)||shape(-x[2],nrow(myds),1);
pmat=exp(myds[,2]#qmat[,1])||exp(myds[,2]#qmat[,1])#qmat[,2]||exp(myds[,2]#qmat[,4])#qmat[,3]||exp(myds[,2]#qmat[,4]);
contrib=pmat[,1]#trds[,1]+pmat[,2]#trds[,2]+pmat[,3]#trds[,3]+pmat[,4]#trds[,4];
loglik=-2#log(sum(contrib#myds[3]));
return(loglik);
finish mylik;
x={0.1,0.8};
opt={0 5 . 3};
call nlpqn(rc,xr,"mylik",x,opt);
quit;
I am here not to ask you about the contruction of pmatrix, qmatrix and selection of contrib, this part might be better of course, but it works and I am happy with it.
What I really need help with is the optimisation part.
In order to mimic the methodology of the R paper I am trying to implement BFGS, with CALL NLPQN opt[4]=3.
Basically I can't manage to get SAS to start the optimisation.
Any help on this part? what am I missing?
below some code to mimic the input dataset:
data mod10; call streaminit(123); whicha=0; do rky=1 to 4; tr1=0; tr2=0; tr3=0; tr4=0; if rky=1 then tr1=1; if rky=2 then tr2=1; if rky=3 then tr3=1; if rky=4 then tr4=1; do rkk=1 to 100; if rand("Uniform")>.5 then do; whicha=whicha+1; timelag=rkk; nocc=1+round(rand("Uniform")*2,1); output; end; end; end; run; proc iml; reset noname; use mod10; read all var {WHICHA TIMELAG NOCC} into myds; read all var {TR1 TR2 TR3 TR4} into trds; close mod10; print myds trds; start mylik(x) global(myds,trds); qmat=shape(-x[1],nrow(myds),1)||shape(x[1],nrow(myds),1)||shape(x[2],nrow(myds),1)||shape(-x[2],nrow(myds),1); pmat=exp(myds[,2]#qmat[,1])||exp(myds[,2]#qmat[,1])#qmat[,2]||exp(myds[,2]#qmat[,4])#qmat[,3]||exp(myds[,2]#qmat[,4]); contrib=pmat[,1]#trds[,1]+pmat[,2]#trds[,2]+pmat[,3]#trds[,3]+pmat[,4]#trds[,4]; loglik=-2#log(sum(contrib#myds[3])); return(loglik); finish mylik; x={0.1,0.8}; opt={0 5 . 3}; call nlpqn(rc,xr,"mylik",x,opt); quit;
I'm on vacation until Tuesday, but here are three hints:
1. Double check your objective function. Usually, we compute SUM(LOG(...)) , not LOG(SUM(...)).
2. Look at this blog post: Two ways to compute maximum likelihood estimates in SAS - The DO Loop
3. Read Ten tips before you run an optimization - The DO Loop (sas.com)
1) I think the problem is you are using EXP() funtion for call nlpqn() . Variable "timelag" has value like >60 is too way big for EXP(). If you make it smaller, you would get convergency result, Like: timelag=rkk; --> timelag=rkk*0.01; 2) And "loglik=-2#log(sum(contrib#myds[3]));" should be "myds[ ,3]" ??
Catch the best of SAS Innovate 2025 — anytime, anywhere. Stream powerful keynotes, real-world demos, and game-changing insights from the world’s leading data and AI minds.