Solved – PROC NLMIXED and PROC LIFEREG not arriving at the same answer for Log-normal survival function

sassurvival

Related question

I have a project where, despite being able to implement some parametric models in LIFEREG, it is somewhat more convenient to do it in NLMIXED. Verifying that this technique works, I tried implementing a pretty bog standard Weibull AFT as follows:

data work.data;
    call streaminit(123);
    do i = 1 to 1000;
        x = rand('BERN', 0.10); *10% exposed;
        t1 = rand('Weibull', 1, 20*(exp(-0.500*x)));
        t2 = t1;
        weight = 1;
        output;
    end;
run;

PROC NLMIXED data=work.data fd;
    parms alpha=1 f0=-1 f1=0;
    bounds alpha>0;
    *Rate of x;
    lam=exp(-(f0*alpha+f1*alpha*x));    

    *Density of x;
    ff1=alpha*lam*t1**(alpha-1)*exp(-lam*t1**alpha);

    *log Likelihood;
    logl=log(ff1);

    *Weighted log Likelihood;
    wlogl=logl*weight;
    model t1~general(wlogl);
    ods exclude iterhistory parameters;
run; quit; run;

proc lifereg data=work.data;    
    model (t1, t2) = x / D=WEIBULL; 
    weight weight; run;

This works swimmingly. Both LIFEREG and NLMIXED produce the same estimates, and all is right in the world. Testing it for a Log-Normal AFT model however using the following modifications doesn't work so well:

PROC NLMIXED data=work.data fd;
    parms sigma=1 g0=-1 g1=0;
    bounds sigma>0;

    *Rate of X;
    mu=exp(g0+g1*x);

    *Density of X;
    fg1 = exp(-0.5*((log(t1)-mu)/sigma)**2)/((t1*(2*CONSTANT('PI'))**0.5)*sigma);

    *log Likelihood;
    logl=log(fg1);

    *Weighted log Likelihood;
    wlogl=logl*weight;
    estimate "RT NDD" exp(g1);
    model x~general(wlogl);
    ods exclude iterhistory parameters;
run;

proc lifereg data=work.data;
    model (t1, t2) = x / D=Lnormal; 
    weight weight;run;

Here, while both models estimate the scale as the same number (1.2676), the estimates for the intercept (g0) and x (g1) are wildly off. In NLMIXED they are 0.8747 and -0.1419 respectively, while in LIFEREG they are 2.3981 and -0.3172.

Any idea what's going on?

Best Answer

Answer:

The rate of x is incorrectly specified as mu=exp(g0+g1*x). Once this is corrected to mu = (g0+g1*x) LIFEREG and NLMIXED arrive at the same answer.

Related Question