froot in SAS/IML

[Salah]

Hello Rick,

I am trying to find a value for the parameter beta that satisfies the equation U_star = 0. I used nlpqn method, but the result returned is the lower limit of the constraint (0.00001). I also tried the FROOT function, but I encountered an error. Additionally, I attempted the NLPLM method, but SAS indicated that I need to provide a vector 2x1. I'm running out of options and would greatly appreciate your help. Thank you.

[ballardw]

Any time you get an error it is a good idea to include the LOG with the code submitted and all of the messages generated. Copy the text from the log, open a text box on the forum using the </> icon and paste the text. The text box prevents the forum from reformatting pasted text and to visually separate the log (or program code) from the discussion text. Also  having the code in the forum post is better than an attachment unless the code is very long.

Best is to include a small data set in the form of working data step code to test the code against as well as to have an idea of the data structures or content that may be causing problems.

[Salah]

Thank you for your reply. The output shows that the estimator based on the Firth method always hits the lower limit of the firth_con; my understanding is that nlpqn finds an optimization, while I want to find a root. This is my primary concern: how to find the solution to the score function=0.

Thank you 

[Rick_SAS · Accepted solution]

Here is an article about the FROOT function, which contains a discussion and example:
A simple way to find the root of a function of one variable - The DO Loop

For your function, here is an example:

proc iml;
      beta = 1;      m = 50;      k = 1;   
      seed  = 2; 

   *-------------------       Firth    -----------------------;
   start Firth_score(y) global(X, m);
      beta = y[1];
      sum_x2 = (X##2)[+];
      /*  Firth  score */
      U_star = - (12*m + 5 - 4*sum_x2/(beta**2)) / (4*beta);
      return(U_star);
   finish;
   
   Call randseed(seed);
   U=j(m,1,0);    X=j(m,1,0); 
   Call randgen(U,"uniform");
   X = Sqrt( (-2*beta**2)# log(1- U) );
   Firth_root = froot("Firth_score", {0.5 1});
   print Firth_root;

Firth_root
0.7115471

You can visualize the Firth score to verify the answer:

   /* visualize to confirm that the Firth_score is zero at root */
   betaVals = T(do(0.5, 1, 0.001)); 
   score = j(nrow(betaVals),1,.); 
   do j = 1 to nrow(betaVals);
     score[j] = Firth_score(betaVals[j]);
   end;
   title "Firth Score as a Function of Beta";
   call series(betaVals, score) grid={x y}  other="refline 0/axis=y;";

[Salah]

Thank you so much!!