# Simplex Model in R

[This article was first published on

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

**Yet Another Blog in Statistical Computing » S+/R**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)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 ### data <- read.dta("/home/liuwensui/Documents/data/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 */

To

**leave a comment**for the author, please follow the link and comment on their blog:**Yet Another Blog in Statistical Computing » S+/R**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

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