why my code doesn't work?

[ddong]

I'm doing a very basic model using SAS/IML for nonlinear optimization by maximizing the likelihood function. Below is my code. It turns out that the model does not converge. The model do converge if I run the it with proc mixed. I was trying to figure out why this happened but I couldn't find out the reason.

proc iml;
use x;
read all into x;
use y;
read all into y;
use nsub;
read all into nsub;
use ntime;
read all into ntime;
use initial;
read all into initial;

nrun=0;
/***loglikelihood funtion***/
start loglik(parameters) global(x,y,nsub,ntime,nrun);
beta=parameters[1:2];
print beta;
rho=parameters[3];
print rho;
sigma2=parameters[4];
print sigma2;

ll=0;
ind=1;
do while (ind<=nsub);
xi=x[(ind-1)*ntime+1:ind*ntime,];
yi=y[(ind-1)*ntime+1:ind*ntime];
pi=constant('PI');

/*remove missing values (monotone missing)*/
if (yi[nrow(yi)]=.) then do; /*if yi is not complete*/
j=nrow(yi);
do while (yi[j-1]=.);
j=j-1;
end;
di=j;
yobs=yi[1:di-1];
end;
else do;
yobs=yi;
end;
ll1=0;
ni=nrow(yobs);
mui=xi[1:ni,]*beta;
help1=0:(ni-1);
help2=shape(help1,ni,ni);
h=help2-t(help2);
h=rho##abs(h); /*elementwise power*/
vi=sigma2*h;
ll1=-(ni/2)*log(2*pi)-0.5*log(det(vi))-0.5*t(yobs-mui)*inv(vi)*(yobs-mui);
/*calculate ll*/
ll=ll+ll1;
print ll;
ind=ind+1;
end;
print ll;
nrun=nrun+1;
file log;
put nrun;
return(ll);
finish loglik;

opt={1 11};
con={. . -1 0 ,
. . 1 . };
call nlpnrr(rc, xr, "loglik", initial, opt, con);
call nlpfdd(maxlik,grad,hessian,"loglik",est);
inf=-hessian;
covar=inv(inf);
var=vecdiag(covar);
stde=sqrt(var);
create result var {est stde};
append;
quit;

[Rick_SAS]

First, if this is a "very basic model" for you, then I'm not sure how much help I can offer, since it doesn't look basic to me!

Second, I have to compliment you. You have clearly done your research.

You are using NLPNRR to match the fact that PROC MIXED uses a ridge-stabilized Newton-Raphson algorithm.

You are using the same ML formula in the PROC MIXED doc.

So what is the source of the problem? I don't see anything obvious, but nonlinear optimization is complicated. The PROC MIXED documentation lists 14 different problems that can occur during the MLE, and many of those are general issues apply to general optimization problems.

You say "that the model does not converge," but you don't give the reason. Does SAS report an error? (If so, what is it? Can you make the error go away by constraining sigma2>1e-8?) Is the gradient not zero at the final estimate (if so, allow more iterations or change the criteria that stop the algorithm). 

One specific suggestion: Take the solution that PROC MIXED gives and plug it into the log-likelihood function that you have define.  Does your IML code give the same gradient and Hessian at that point as does PROC MIXED? If not, your function and PROC MIXED are computing different things.  (The LL is only defined up to a constant, so the LL might not match.)

A related suggestion: if you use the solution from PROC MIXED as an initial condition, does the IML code converge?  If so, maybe the intial condition is bad.

Rick

PS. One small performance suggestion: You don't need a loop to eliminate the missing values. Use the LOC function to find the nonmissing values, then extract them. For example, I think your IF-THEN...DO WHILE code is equivalent to yobs = yi[loc(y^=.)]; For more on removing missing values, see http://blogs.sas.com/content/iml/2010/09/15/removing-observations-with-missing-values/

[ddong]

Hi Rick,

Thank you very much for your reply. It is really appreciated.

Sas did not report any error when I ran the code, it just kept runing the interation for the optimization but never stopped. And I did use the solution from Proc Mixed as the intial condition, but seems it didn't help. That's where I'm really confused.

BTW, thanks for your suggestion of dealing with the missing values.

Dongli

[ddong]

BTW, i was trying to monitor the changes of the estimates so I specified print beta... within the loglik function. It seems that the estimates did not update from each interation.

[Rick_SAS]

Perhaps the likelihood function is extremely flat and so each new step is only 1e-12 away from the old one. A situation like that would take a long time to converge.

If that's the case, does it really matter whether you have the exact optimum or are just at a parameter value that is nearly optimum?  You might relax the convergence criterion so that the algorithm stops when it reaches the approximate optimum.

[ddong]

Thanks Rick. This is really helpful.

[Logie]

Hi Rick

I wonder if you can help. What could be causing SAS Enterprise Guide to return the following error when using the Query Builder tab?

Many thanks

Gwin

****************************************************************************************************************************************

System.Xml.XmlException

Root element is missing.

------------------------------Technical Information Follows ------------------------------

System.Xml.XmlException: Rootelement is missing.

   atSystem.Xml.XmlTextReaderImpl.Throw(Exception e)

   atSystem.Xml.XmlTextReaderImpl.ThrowWithoutLineInfo(String res)

   atSystem.Xml.XmlTextReaderImpl.ParseDocumentContent()

   atSystem.Xml.XmlTextReaderImpl.Read()

   atSystem.Xml.XmlLoader.Load(XmlDocument doc, XmlReader reader, BooleanpreserveWhitespace)

   atSystem.Xml.XmlDocument.Load(XmlReader reader)

   atSAS.SharedUI.Controls.ExpressionBuilderControl.ReadInFile(String strPath)

   atSAS.SharedUI.Controls.ExpressionBuilderControl.AddFunctions(String filter)

   atSAS.EG.DataAccess.EGExpressionBuilder.EGExpressionControl.AddFunctions(Stringfilter)

   atSAS.EG.ProjectElements.Forms.FilterAdvancedPage..ctor(QueryBuilderQueryBuilder, ShowComputedColumns ComputedColumns)

   atSAS.EG.ProjectElements.Forms.FilterWizard.FilterDialog(FilterTreeEnumfiltertype)

   atSAS.EG.ProjectElements.Forms.QueryBuilder.AddFilter(FilterTree tree, FilterNodeparent)

   atSAS.EG.ProjectElements.Forms.QueryBuilder.OnNewFilterRequested(Object sender,FilterListEventArgs e)

   atSAS.EG.ProjectElements.Forms.QueryBuilder.FilterList_FilterEvent(Object sender,FilterListEventArgs e)

   atSAS.SharedUI.QueryFilters.FilterList.cmdNewFilter_Click(Object sender,EventArgs e)

   atSystem.Windows.Forms.Control.OnClick(EventArgs e)

   atSystem.Windows.Forms.Button.OnClick(EventArgs e)

   atSystem.Windows.Forms.Button.OnMouseUp(MouseEventArgs mevent)

   atSAS.SharedUI.Controls.SASImagePushButton.OnMouseUp(MouseEventArgs e)

   atSystem.Windows.Forms.Control.WmMouseUp(Message& m, MouseButtons button,Int32 clicks)

   at System.Windows.Forms.Control.WndProc(Message&m)

   atSystem.Windows.Forms.ButtonBase.WndProc(Message& m)

   atSystem.Windows.Forms.Button.WndProc(Message& m)

   atSAS.SharedUI.Controls.SASImagePushButton.WndProc(Message& m)

   atSystem.Windows.Forms.Control.ControlNativeWindow.OnMessage(Message& m)

   atSystem.Windows.Forms.Control.ControlNativeWindow.WndProc(Message& m)

   atSystem.Windows.Forms.NativeWindow.Callback(IntPtr hWnd, Int32 msg, IntPtrwparam, IntPtr lparam)

[Rick_SAS]

Sorry, but I have no idea. Try posting your question at

http://communities.sas.com/community/sas_enterprise_guide