PRJMultinomial—predicted probs for a multinomial logit
Posted: Thu Jul 10, 2008 2:01 pm
This computes the predicted probabilities and marginal effects for the choices in a multinomial logit model given an input set of explanatory variables.
Code: Select all
*
* @PRJMultinomial( options )
* computes predicted probabilities and marginal effects for choices in a
* multinomial logit. This should be only be used after doing DDV with
* TYPE=MULTINOMIAL.
*
* Options:
* XVECTOR=vector of test values for the explanatory variables. Make sure you include
* 1 for the CONSTANT.
* LABELS=VECTOR[STRINGS] with labels for the choices. These should correspond to the
* numeric codings from smallest to largest.
* TITLE=title for output
*
* Revision Schedule:
* 07/2008 Written by Tom Doan, Estima.
*
procedure prjmultinomial
option vector xvector
option vect[strings] labels
option string title
*
local real p
local vect ez pz dg m
local rect betamat
local integer i j row
*
if %rows(%betasys)==0.or.%clock(%rows(%betasys),%nreg)<>%nreg {
disp "@PRJMULTINOMIAL must follow DDV(TYPE=MULTINOMIAL)"
return
}
*
if .not.%defined(xvector) {
disp "Syntax: @PRJMULTINOMIAL(XVECTOR=xvector,other options)"
return
}
*
* Reformat %betasys into an %NREG x (choices - 1) array
*
compute betamat=%vectorect(%betasys,%nreg)
*
* Compute the probabilities of the choices (other than 1)
*
compute ez=%exp(tr(xvector)*betamat)
compute pz=ez/(1+%sum(ez))
compute dg=betamat*pz
report(action=define)
compute row=0
if %defined(title)
report(atrow=(row=row+1),atcol=1,span) title
do i=1,%rows(ez)+1
*
* Handle the first choice by taking the remainders
*
if i==1 {
compute p=1-%sum(pz)
compute m=-p*dg
}
else {
compute p=pz(i-1)
compute m=p*(%xcol(betamat,i-1)-dg)
}
compute row=row+1
if %defined(labels)
report(atrow=row,atcol=1) "Choice: "+labels(i)
else
report(atrow=row,atcol=1) "Choice "+i
report(atrow=(row=row+1),atcol=1) "Predicted P" p
do j=1,%nreg
if %eqnreglabels(0)(j)=="Constant"
next
report(atrow=(row=row+1),atcol=1) "Marginal "+%eqnreglabels(0)(j) m(j)
end do j
compute row=row+1
end do i
report(action=show)
end