Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

R CODE

```library(simplexreg)
library(foreign)

### http://fmwww.bc.edu/repec/bocode/k/k401.dta ###

mdl <- simplexreg(prate ~ mrate + totemp + age + sole|mrate + totemp + age + sole, type = "hetero", link = "logit", data = data, subset = prate < 1)

summary(mdl)
```

R OUTPUT

```simplexreg(formula = prate ~ mrate + totemp + age + sole | mrate + totemp +
age + sole, data = data, subset = prate < 1, type = "hetero", link = "logit")

standard Pearson residuals:
Min      1Q  Median      3Q     Max
-6.1724 -0.5369  0.0681  0.5379  2.2987

Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)  8.817e-01  4.036e-02  21.848  < 2e-16 ***
mrate        2.710e-01  4.880e-02   5.553 2.81e-08 ***
totemp      -8.454e-06  1.164e-06  -7.266 3.70e-13 ***
age          2.762e-02  2.702e-03  10.225  < 2e-16 ***
sole         1.079e-01  4.684e-02   2.304   0.0212 *

Coefficients (dispersion model with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.668e+00  5.395e-02  30.918  < 2e-16 ***
mrate       8.775e-01  4.472e-02  19.621  < 2e-16 ***
totemp      7.432e-06  1.434e-06   5.182  2.2e-07 ***
age         2.816e-02  3.242e-03   8.688  < 2e-16 ***
sole        7.744e-01  5.966e-02  12.981  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-likelihood:  2370,  p-value: 0.4693177
Deviance: 2711
Number of Fisher Scoring iterations:  20
```

SAS CODE & OUTPUT FOR COMPARISON

```proc nlmixed data = one tech = trureg maxiter = 100;
parms b0 = 0  b1 = 0  b2 = 0  b3 = 0  b4 = 0
c0 = 2  c1 = 0  c2 = 0  c3 = 0  c4 = 0 ;
xb = b0 + b1 * mrate + b2 * totemp + b3 * age + b4 * sole;
xc = c0 + c1 * mrate + c2 * totemp + c3 * age + c4 * sole;
mu_xb = 1 / (1 + exp(-xb));
s2 = exp(xc);
v = (prate * (1 - prate)) ** 3;
d = (prate - mu_xb) ** 2 / (prate * (1 - prate) * mu_xb ** 2 * (1 - mu_xb) ** 2);
lh = (2 * constant('pi') * s2 * v) ** (-0.5) * exp(-(2 * s2) ** (-1) * d);
ll = log(lh);
model prate ~ general(ll);
run;
/*
Standard
Parameter   Estimate      Error     DF   t Value   Pr > |t|    Alpha
b0            0.8817    0.03843   2711     22.94     <.0001     0.05
b1            0.2710    0.04540   2711      5.97     <.0001     0.05
b2          -8.45E-6    1.35E-6   2711     -6.26     <.0001     0.05
b3           0.02762   0.002588   2711     10.67     <.0001     0.05
b4            0.1079    0.04792   2711      2.25     0.0244     0.05
c0            1.6680    0.05490   2711     30.38     <.0001     0.05
c1            0.8775    0.07370   2711     11.91     <.0001     0.05
c2          7.431E-6   1.935E-6   2711      3.84     0.0001     0.05
c3           0.02816   0.003224   2711      8.73     <.0001     0.05
c4            0.7744    0.06194   2711     12.50     <.0001     0.05
*/
```  