Good Day.
I'm trying to determine how SAS calculates the predicted values in proc lifereg. The calculation is described here: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_lifereg_sec... . The default value of p is .5. The code below shows an example.
data work.SomeStuff (drop = i);
call streaminit(666);
do i = 1 to 100000;
PredictorVariable = rand('uniform');
LinearPartOfTheMean = 2 * PredictorVariable + 3;
/* Create a gamma random variable */
GammaVariable = rand('gamma', exp(LinearPartOfTheMean));
output;
end;
run;
proc lifereg data = work.SomeStuff
outest = work.GammaEstimates (keep = Intercept PredictorVariable _scale_ _shape1_
rename = (Intercept = GammaIntercept PredictorVariable = GammaPredictorVariable _scale_ = GammaScale _shape1_ = GammaShape));
model GammaVariable = PredictorVariable / distribution = gamma;
output out = work.SomeStuff02 (drop = _prob_) predicted = GammaPrediction;
run;
data work.SomeStuff02 (drop = GammaIntercept GammaPredictorVariable GammaScale GammaShape);
set work.SomeStuff02;
if _n_ = 1 then set work.GammaEstimates;
GammaPrediction02 = exp(GammaIntercept + GammaPredictorVariable * PredictorVariable);
Ratio = GammaPrediction / GammaPrediction02;
run;
proc print data = work.SomeStuff02 (obs = 10);
where GammaPrediction < GammaPrediction02;
run;
proc means data = work.SomeStuff02 min max;
var Ratio;
run;
proc print data = work.GammaEstimates;
run;
The SAS documentation indicates that I need to add the scale parameter multiplied by the value of u_p, which for the gamma distribution is computed numerically, to what is exponentiated in GammaPrediction02. How is the value of u_p computed? Moreover, since the scale parameter in the example is positive, the value of u_p must be negative to achieve a value less than 1 for the value of Ratio. However, this is impossible since u_p is a quantile from a gamma distribution which has only positive support.
Thank you for your time and attention in this matter.
