# Yet Another Blog in Statistical Computing

I can calculate the motion of heavenly bodies but not the madness of people. -Isaac Newton

## Two-Part Fractional Model

1. Two-step estimation method:
– SAS code

```data one;
set data.sme;
_id_ + 1;
if y = 0 then y2 = 0;
else y2 = 1;
run;

proc logistic data = one desc;
model y2 = x1 - x7;
run;

proc nlmixed data = one tech = congra;
where y2 = 1;
parms b0 =  2.14 b1 = -0.01 b2 = -0.55 b3 = -0.14 b4 = -3.47
b5 =  0.01 b6 = -0.01 b7 = -0.19 s = 0.1;
_xb_ = b0 + b1 * x1 + b2 * x2 + b3 * x3 + b4 * x4 +
b5 * x5 + b6 * x6 + b7 * x7;
_mu_ = 1 / (1 + exp(-_xb_));
lh = pdf('normal', y, _mu_, s);
ll = log(lh);
model y ~ general(ll);
predict ll out = out1(keep = _id_ pred rename = (pred = ln_like1));
run;
```

– Outputs

```*** output for logit component ***
Model Fit Statistics
Intercept
Intercept            and
Criterion          Only     Covariates
AIC            4997.648       4066.398
SC             5004.043       4117.551
-2 Log L       4995.648       4050.398

Analysis of Maximum Likelihood Estimates
Standard          Wald
Parameter    DF    Estimate       Error    Chi-Square    Pr > ChiSq
Intercept     1     -9.3586      0.4341      464.8713        <.0001
x1            1     -0.0595      0.0341        3.0445        0.0810
x2            1      1.7644      0.1836       92.3757        <.0001
x3            1      0.5994      0.0302      392.7404        <.0001
x4            1     -2.5496      0.5084       25.1488        <.0001
x5            1    -0.00071     0.00128        0.3094        0.5780
x6            1    -0.00109     0.00267        0.1659        0.6838
x7            1     -1.6359      0.2443       44.8427        <.0001

*** output for the positive component ***
Fit Statistics
-2 Log Likelihood                 -264.5
AIC (smaller is better)           -246.5
AICC (smaller is better)          -246.3
BIC (smaller is better)           -201.4

Parameter Estimates
Standard
Parameter   Estimate      Error     DF   t Value   Pr > |t|
b0            1.7947     0.3401   1116      5.28     <.0001
b1          -0.01216    0.02737   1116     -0.44     0.6568
b2           -0.4306     0.1437   1116     -3.00     0.0028
b3           -0.1157    0.02287   1116     -5.06     <.0001
b4           -2.9267     0.4656   1116     -6.29     <.0001
b5          0.004092   0.001074   1116      3.81     0.0001
b6          -0.00838   0.002110   1116     -3.97     <.0001
b7           -0.1697     0.1977   1116     -0.86     0.3909
s             0.2149   0.004549   1116     47.24     <.0001
```

2. Joint estimation method:
– SAS code

```proc nlmixed data = one tech = congra maxiter = 1000;
parms b10 = -9.3586 b11 = -0.0595 b12 =  1.7644 b13 =  0.5994 b14 = -2.5496
b15 = -0.0007 b16 = -0.0011 b17 = -1.6359
b20 =  0.3401 b21 =  0.0274 b22 =  0.1437 b23 =  0.0229 b24 =  0.4656
b25 =  0.0011 b26 =  0.0021 b27 =  0.1977  s  =  0.2149;
logit_xb = b10 + b11 * x1 + b12 * x2 + b13 * x3 + b14 * x4 +
b15 * x5 + b16 * x6 + b17 * x7;
nls_xb = b20 + b21 * x1 + b22 * x2 + b23 * x3 + b24 * x4 +
b25 * x5 + b26 * x6 + b27 * x7;
p1 = 1 / (1 + exp(-logit_xb));
if y = 0 then ll = log(1 - p1);
else ll = log(p1) + log(pdf('normal', y, 1 / (1 + exp(-nls_xb)), s));
model y ~ general(ll);
run;
```

– Outputs

```             Fit Statistics
-2 Log Likelihood                 3785.9
AIC (smaller is better)           3819.9
AICC (smaller is better)          3820.0
BIC (smaller is better)           3928.6

Parameter Estimates
Standard
Parameter   Estimate      Error     DF   t Value   Pr > |t|
b10          -9.3586     0.4341   4421    -21.56     <.0001
b11         -0.05948    0.03408   4421     -1.75     0.0810
b12           1.7645     0.1836   4421      9.61     <.0001
b13           0.5994    0.03024   4421     19.82     <.0001
b14          -2.5496     0.5084   4421     -5.01     <.0001
b15         -0.00071   0.001276   4421     -0.56     0.5784
b16         -0.00109   0.002673   4421     -0.41     0.6836
b17          -1.6359     0.2443   4421     -6.70     <.0001
b20           1.7633     0.3398   4421      5.19     <.0001
b21         -0.01115    0.02730   4421     -0.41     0.6830
b22          -0.4363     0.1437   4421     -3.04     0.0024
b23          -0.1139    0.02285   4421     -4.98     <.0001
b24          -2.8755     0.4643   4421     -6.19     <.0001
b25         0.004091   0.001074   4421      3.81     0.0001
b26         -0.00839   0.002109   4421     -3.98     <.0001
b27          -0.1666     0.1975   4421     -0.84     0.3991
s             0.2149   0.004550   4421     47.24     <.0001
```

As shown above, both estimation methods give similar parameter estimates. The summation of log likelihood from the 2-step estimation models is exactly equal to the log likelihood from the joint estimation model.