Convergence Problems with proc lifereg

[Salah]

Hello 

I am running long simulations and using proc lifereg.  I need to compare the results based on different censoring schemes and different sample size m. One of the problems that I faced is that the Mean squared error (MSE) is not decreasing as sample size increases which is wiered. I thought that maybe I have a convergence problem. So I searched  that and the solution that I found on the web is  is to fit a log logistic distribution, then use the resulting parameter estimates as  initial values in a subsequent fit of a model with the Weibull distribution.

My problem is that I have simulation and in each simulation I get a different data set therefore different parameter estimates using log logistic. How would I solve this problem. 

Thank you for your time and help.

[Rick_SAS]

Let me answer the MSE question in terms of least squares regression, which is easier to describe and more familiar. The standard errors of the parameter estimates should get smaller as the sample size increases.  (For linear regression, the standard errors are sqrt(MSE/n).) However, the RMSE = sqrt(MSE) is a measure of the noise component that you add to the linear predictor. It is fixed and does not get smaller with larger sample size. 

> fit a log logistic distribution, then use the resulting parameter estimates as initial values 

OK, I am not familiar with this method, but I understand the main idea of fitting a reduced or associated model.

> My problem is that I have simulation and in each simulation I get a different data set

> therefore different parameter estimates using log logistic. How would I solve this problem?

I'm not sure what problem you have. Yes, each simulation creates a different sample and therefore different LL param estimates. Use those estimates for the LIFEREG analysis of that sample. 

[Salah]

Thank you Rick

>I'm not sure what problem you have. Yes, each simulation creates a different sample and therefore different LL param estimates. Use those estimates for the LIFEREG analysis of that sample. 

The problem is that I can't stop the simulation to copy the results from the first proc lifprog using log-logistic  and enter it in the second proc lifprog using weibull. Is there away I can save the results of 

proc lifereg data=AFT;
model T1*U_2(0) = Z1/ distribution = llogistic;
run;

Let us call them intercept1 , initial1, scale1 then use them below like this

proc lifereg data=AFT outest=AFT_out covout ;
model T1*U_2(0) = Z1 / itprint dist=WEIBULL
intercept=intercept1 initial= initial1   scale=scale1;
run;

I hope what I am saying makes sense

Thank you.

[IanWakeling]

I think the trick you are looking for is to use PROC SQL to get the parameter estimates into macro variables.  The estimates themselves can be read from the OUTEST data set in LIFEREG.

proc lifereg data=AFT outest=est;
  model T1*U_2(0) = Z1 / itprint distribution = llogistic;
run;

proc sql noprint;
   select  intercept format=best20.,
           Z1        format=best20.,
           _SCALE_   format=best20.  into
           :p1 TRIMMED, :p2 TRIMMED, :p3 TRIMMED from est;
quit;

%put &=p1 &=p2 &=p3;

then in the next model you can specify

intercept=&p1  initial=&p2  scale=&p3

as options on the model statement.

[Salah]

Thank you so much Ian

This is exactly what I wont. However I received an error message which I couldn't figure it out or know how to fix it. Please see the attached file.

Thank you

/**AFT wiebul  using Progressive TypeII censoring **/

goptions reset=all;
%macro ODSOff(); /* Call prior to BY-group processing */
ods graphics off;
ods exclude all;
ods noresults;
%mend;

%macro ODSOn(); /* Call after BY-group processing */
ods graphics on;
ods exclude none;
ods results;
%mend;

%odsoff;
proc iml;

k=100; seed=0;    

B1=0;

scheme=16;n=300; m=280; R=J(m,1,0); do i=1 to m by 14; R[i]=1; end;

******************************;


odds=exp(B1);
lambda=0.5;B0=1;


**** Defining the variables for the Inference part ****;

Beta1=J(k,1,0);

Bias_B1=J(K,1,0);

MSE_B1=J(K,1,0);Var_B1=J(K,1,0); 

sq1=J(K,1,0);    Lower_bound=J(K,1,0);    Upper_bound=J(K,1,0);    Test1=J(K,1,0);   odd1=J(K,1,0);    
MSE_ods=J(K,1,0); Bias_ods=J(K,1,0);     


**********************************************************;

Do N1=1 to K; *simulation starts here*;


** Simulation of progressive censoring ****;

hu=J(n,1,0);

Ww=J(m,1,0);X=J(m,1,0);UU1=J(m,1,0);U2=J(m,1,0);T1=J(n,1,0); V=J(n,1,0);
Y1=J(m,1,0);Z1=J(n,1,0);Z2=J(n,1,0);

Ee=J(m,1,0);Vv=J(m,1,0);ll=J(m,1,0);
 

** Creation of the progressive data from Uniform;
************************************************;	

do i=1 to m;

						
				Ww[i]=Uniform(seed);
				ll[i]=i;

						sum=0;	
						do j=m-i+1 to m;
						sum=(R[j]+sum);
						end;

				Ee[i]=1/(i+sum);
				Vv[i]=Ww[i]**Ee[i];
				
end;

	do i=1 to m;

				prod=1;
				endloop=m-i+1 ; 
						do k1=m to endloop by -1;
						
							prod=Vv[k1]*prod;
							
						end;

				UU1[i]=1-(prod);** the U's are the required progressively type II
	right censored sample from U(0,1);

	end;

	U_obs=UU1||J(m,1,1);** progressive censored data of size m from U(0,1) with idex=1 to indicate the event occurs;



** I need the data to have a size n, so I need to generate more data from Uniform(0,1) with size of n-m; 

UU0=J(n-m,1,0);

do j=1 to n-m;

	UU0[j]=Uniform(seed);

end;

U_Cen=UU0||J(n-m,1,0); * These are the censored data with idex 0;



U=J(n,2,0);

U=U_obs//U_cen;


U_1=U[,1];* This is the first column of the matrix U which represents the data from Uniform both(observed + censored); 
U_2=U[,2];* This is the index (1:failure  OR  0:Censored);



** Simulate SRS data from Weibull **;

do jj=1 to n;

hu[jj]=Uniform(1111);

end;

Z1=quantile('Normal', hu);


C=B0+B1#Z1;


T1=exp(C)#(-log(U_1))##lambda;	/*Weibull mixed between progressive and ordinary data**;



/* (1) To use the values of my r.v. in proc lifereg, we need to write them to a SAS data set using CREATE and APPEND 
statements. These statements create data set AFT;*/


create AFT var { "U_obs", "U_Cen", "U" , "U_1" ,"U_2","T1", "Z1"};
append;
close AFT;

* (2) Next, the lifereg procedure is called from within IML by using SUBMIT and ENDSUBMIT statements. ;


submit;
proc lifereg data=AFT;
   model T1*U_2(0) = Z1/ distribution = llogistic;
run;
proc sql noprint;
   select  intercept format=best20.,
           Z1        format=best20.,
           _SCALE_   format=best20.  into
           :p1 TRIMMED, :p2 TRIMMED, :p3 TRIMMED from est;
quit;

%put &=p1 &=p2 &=p3;

proc lifereg data=AFT outest=AFT_out covout ;
model T1*U_2(0) = Z1 / itprint dist=WEIBULL 
                       intercept=&p1  initial=&p2  scale=&p3;
run;

data new_AFT;
set AFT_out;
keep  Z1;
run;

proc transpose data=new_AFT out=idnumber name=Test
     prefix=sn;
run;
endsubmit;


* (3) To read the results back to SAS/IML we write "use", "read", and "close";

Use idnumber;
read all var {Sn1 Sn2 Sn3 Sn4};
Close idnumber;

Beta1[N1]=Sn1;

Var_B1[N1]=Sn3;

Use AFT;
read all var {U_2};
Close AFT;


********** Inference ***********;

Beta1[N1]=Sn1;

odd1[N1]=exp(Sn1);

Var_B1[N1]=Sn3;

MSE_B1[N1]=(Beta1[N1]-B1)**2;

MSE_ods[N1]=(odd1[N1]-odds)**2;

end; *end of simulation loop (NN);

%odson

** Estimation of Beta1 ** ;
************************* ;

Beta1_hat=(Beta1[:]);	
Bias_B1=abs(Beta1[:]-B1);
Bias_Beta1=round(Bias_B1,0.00001);

MSE_Beta1=round(MSE_B1[:],0.00001);

** Estimation of hazards ration OR Odds ** ;
****************************************** ;

odd1_hat=(odd1[:]);
Biasodd1=abs(odd1_hat-Odds);
Bias_odd1=round(Biasodd1,0.00001);

MSE_Odd1=round(MSE_ods[:],0.00001);




print Bias_Beta1  MSE_Beta1   Bias_odd1  MSE_Odd1;
 
  
 
  Quit;

[Rick_SAS · Accepted solution]

You didn't copy Ian's code correctly:

proc lifereg data=AFT outest=est;

[Salah]

You are right Rick

Thank you so much for your swift reply.

[Salah]

Thank you Ian

I am so grateful for your trick. It helps me more than you can imagine. AMAZING trick!!

[Rick_SAS]

In addition to the question you asked, you might want to read the article, "Eight tips to make your simulation run faster."

The first tip is to vectorize your IML code by replacing loops with vector operations. For example, in your program, here are two easy loops to replace.

/*
      sum=0;
      do j=m-i+1 to m;
         sum=(R[j]+sum);
      end;
*/
      sum = sum(R[(m-i+1):m]); 

...
/*
      prod=1;
      endloop=m-i+1;
      do k1=m to endloop by -1;
         prod=Vv[k1]*prod;
      end;
*/
      prod = prod(Vv[m:(m-i+1)]);

[Salah]

Yes sure sure thank you so much

[Rick_SAS]

You might want to mark Ian's answer as the correct solution.