(
newcommand{xb}{{bf x}}
newcommand{betab}{boldsymbol{beta}})I talk about mypoisson1, which computes Poisson-regression ends in Mata. The code in mypoisson1.ado is remarkably much like the code in myregress11.ado, which computes atypical least-squares (OLS) ends in Mata, as I mentioned in Programming an estimation command in Stata: An OLS command utilizing Mata.
I construct on earlier posts. I exploit the construction of Stata applications that use Mata work capabilities that I mentioned beforehand in Programming an estimation command in Stata: A primary ado-command utilizing Mata and Programming an estimation command in Stata: An OLS command utilizing Mata. Try to be aware of Poisson regression and utilizing optimize(), which I mentioned in Programming an estimation command in Stata: Utilizing optimize() to estimate Poisson parameters.
That is the nineteenth submit within the collection Programming an estimation command in Stata. I like to recommend that you simply begin firstly. See Programming an estimation command in Stata: A map to posted entries for a map to all of the posts on this collection.
A poisson command with Mata computations
The Stata command mypoisson1 computes the ends in Mata. The syntax of the mypoisson1 command is
mypoisson1 depvar indepvars [if] [in] [, noconstant]
the place indepvars can comprise time-series variables. mypoisson1 doesn’t permit for issue variables as a result of they complicate this system. I talk about these problems, and current options, in my subsequent submit.
Within the the rest of this submit, I talk about the code for mypoisson1.ado. I like to recommend that you simply click on on the file identify to obtain the code. To keep away from scrolling, view the code within the do-file editor, or your favourite textual content editor, to see the road numbers.
Code block 1: mypoisson1.ado
*! model 1.0.0 31Jan2016
program outline mypoisson1, eclass sortpreserve
model 14.1
syntax varlist(numeric ts min=2) [if] [in] [, noCONStant ]
marksample touse
gettoken depvar indepvars : varlist
_rmcoll `indepvars', `fixed' forcedrop
native indepvars "`r(varlist)'"
tempname b V N rank
mata: mywork("`depvar'", "`indepvars'", "`touse'", "`fixed'", ///
"`b'", "`V'", "`N'", "`rank'")
if "`fixed'" == "" {
native cnames "`indepvars' _cons"
}
else {
native cnames "`indepvars'"
}
matrix colnames `b' = `cnames'
matrix colnames `V' = `cnames'
matrix rownames `V' = `cnames'
ereturn submit `b' `V', esample(`touse') buildfvinfo
ereturn scalar N = `N'
ereturn scalar rank = `rank'
ereturn native cmd "mypoisson1"
ereturn show
finish
mata:
void mywork( string scalar depvar, string scalar indepvars,
string scalar touse, string scalar fixed,
string scalar bname, string scalar Vname,
string scalar nname, string scalar rname)
{
actual vector y, b, mo_v, cv
actual matrix X, V, Cm
actual scalar n, p, rank, ko
y = st_data(., depvar, touse)
n = rows(y)
X = st_data(., indepvars, touse)
if (fixed == "") {
X = X,J(n, 1, 1)
}
p = cols(X)
S = optimize_init()
optimize_init_argument(S, 1, y)
optimize_init_argument(S, 2, X)
optimize_init_evaluator(S, &plleval2())
optimize_init_params(S, J(1, p, .01))
b = optimize(S)
V = optimize_result_V_oim(S)
rank = p - diag0cnt(invsym(V))
st_matrix(bname, b)
st_matrix(Vname, V)
st_numscalar(nname, n)
st_numscalar(rname, rank)
}
void plleval2(actual scalar todo, actual vector b, ///
actual vector y, actual matrix X, ///
val, grad, hess)
{
actual vector xb
xb = X*b'
val = sum(-exp(xb) + y:*xb - lnfactorial(y))
}
finish
mypoisson1.ado has the construction of Stata applications that compute their ends in Mata, which I mentioned in Programming an estimation command in Stata: A primary ado-command utilizing Mata and Programming an estimation command in Stata: An OLS command utilizing Mata. Strains 1–35 outline the ado-command mypoisson1. Strains 37–83 outline the Mata work operate mywork() utilized in mypoisson1 and the evaluator operate plleval2() utilized in mywork().
The ado-command mypoisson1 has 4 components:
- Strains 5–13 parse what the consumer typed, establish the pattern, drop collinear variables from the record of unbiased variables, and create short-term names for Stata objects returned by our Mata work operate.
- Strains 15–16 name the Mata work operate.
- Strains 18–31 submit the outcomes returned by the Mata work operate to e().
- Line 33 shows the outcomes.
The Mata work operate mywork() has 4 components.
- Strains 39–42 parse the arguments.
- Strains 45–47 declare vectors, matrices, and scalars which can be native to mywork().
- Strains 49–65 compute the outcomes.
- Strains 67–70 copy the computed outcomes to Stata, utilizing the names that had been handed in arguments.
Now, I talk about the ado-code in some element. Strains 2–35 are nearly the identical as strains 2–35 of myregress11, which I mentioned in Programming an estimation command in Stata: An OLS command utilizing Mata. That myregress11 handles issue variables however mypoisson1 doesn’t causes a lot of the variations. That myregress11 shops the residual levels of freedom however mypoisson1 doesn’t causes two minor variations.
myregress11 makes use of the matrix inverter to deal with circumstances of collinear unbiased variables. Within the Poisson-regression case mentioned right here, collinear unbiased variables outline an unconstrained drawback and not using a distinctive answer, so optimize() can not converge. mypoisson1 drops collinear unbiased variables to keep away from this drawback. I talk about a greater answer in my subsequent submit.
Strains 10–11 use _rmcoll to drop the collinear unbiased variables and retailer the record of linearly unbiased variables within the native macro indepvars.
Now, I talk about the Mata work operate mywork() in some element. The mywork() outlined on strains 39–71 is remarkably much like the mywork() operate outlined on strains 39–72 of myregress11.ado, mentioned in Programming an estimation command in Stata: An OLS command utilizing Mata. In mypoisson1.ado, strains 57–65 use optimize() to compute the outcomes. In myregress11.ado, strains 58–64 use matrix computations to compute the outcomes.
Strains 73–81 outline the evaluator operate plleval2(), which is utilized by optimize() to compute the outcomes, as I mentioned in Programming an estimation command in Stata: Utilizing optimize() to estimate Poisson parameters.
Examples 1 and a couple of illustrate that mypoisson1 produces the identical outcomes as poisson.
Instance 1: mypoisson1
(Makes use of accident3.dta)
. clear all
. use accident3
. mypoisson1 accidents cvalue youngsters site visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66874
Iteration 2: f(p) = -555.81708
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943554 -.5174188
youngsters | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
site visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082077
_cons | .5743543 .2839519 2.02 0.043 .0178187 1.13089
------------------------------------------------------------------------------
Instance 2: poisson
poisson accidents cvalue youngsters site visitors
Iteration 0: log probability = -555.86605
Iteration 1: log probability = -555.8154
Iteration 2: log probability = -555.81538
Poisson regression Variety of obs = 505
LR chi2(3) = 340.20
Prob > chi2 = 0.0000
Log probability = -555.81538 Pseudo R2 = 0.2343
------------------------------------------------------------------------------
accidents | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
youngsters | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506594
site visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .574354 .2839515 2.02 0.043 .0178193 1.130889
------------------------------------------------------------------------------
Completed and undone
I mentioned mypoisson1, which computes Poisson-regression ends in Mata. I highlighted how related the code in mypoisson1.ado is to the code in myregress11.ado.
I additionally mentioned that mypoisson1 drops collinear unbiased variables and famous that mypoisson1 doesn’t deal with issue variables. In my subsequent submit, I talk about an answer to the issues brought on by collinear unbiased variables and talk about a command that handles issue variables.
