[R] Help Computing Probit Marginal Effects

[Peter Ehlers]

On 2010-02-27 1:38, (Ted Harding) wrote:
> On 27-Feb-10 03:52:19, Cardinals_Fan wrote:
>>
>> Hi,  I am a stata user trying to transition to R.  Typically I
>> compute marginal effects plots for (example) probit models by
>> drawing simulated betas by using the coefficient/standard error
>> estimates after I run a probit model. I then use these simulated
>> betas to compute first difference marginal effects.  My question
>> is, can I do this in R?  Specifically, I was wondering if anyone
>> knows how R stores the coefficient/standard error estimates after
>> you estimate the model?  I assume it's  vector, but what is it
>> called?
>>
>> Cheers
>> --
>
> Here is an example which sets up (X,Y) data using a probit mechaism,
> then fits a probit model, and then extracts the information which
> you seek.
>
>    set.seed(54321)
>    X<- 0.2*(-10:10)
>    U<- rnorm(21)
>    Y<- 1*(U<= X)  ## binary outcome 0/1, = 1 if N(0,1)<= X
>    GLM<- glm(Y ~ X, family=binomial(link=probit)) ## fit a probit
>    Coef<- summary(GLM)$coef  ## apply summary() to the fit
>
> GLM is a list with a large number of components: enter the command
>
>    str(GLM)
>
> and have a look at what you get! Only a few of these are displayed
> when you apply print() to it:
>
>    print(GLM)
>    # Call:  glm(formula = Y ~ X, family = binomial(link = probit))
>    # Coefficients:
>    # (Intercept)            X
>    #     0.08237      0.56982
>    #
>    # Degrees of Freedom: 20 Total (i.e. Null);  19 Residual
>    # Null Deviance:      29.06
>    # Residual Deviance: 23.93        AIC: 27.93
>
> Note that you do *not* get Standard Errors from this.
>
> However, all the information in GLM is available for processing
> by other functions. In particular, summary(GLM) produces another
> list with several components -- have a look at the output from
>
>    str(summary(GLM))
>
> One of these components (listed near the end of this output)
> is "coef", and it can be accessed as summary(GLM)$coef as in the
> above command
>
>    Coef<- summary(GLM)$coef
>
> This is a matrix (in this case 2 named rows, 4 named columns):
>
>    Coef
>    #              Estimate Std. Error   z value   Pr(>|z|)
>    # (Intercept) 0.0823684  0.2974595 0.2769063 0.78185207
>    # X           0.5698200  0.2638657 2.1595076 0.03081081
>
> So there is one row for each coefficient in the model (here 2,
> one for Intercept, one for variable X), and four columns
> (for the Estimate itself of the coefficient, for its Standard
> Error, for the z-value (Est/SE), and for the P-value).
>
> Hence you can access the estimates as
>
>    Coef[,1]   # (the first column of the matrix)
>    # (Intercept)           X
>    #   0.0823684   0.5698200
>
> and their respective SEs as
>
>    Coef[,2]   # (the second column of the matrix)
>    # (Intercept)           X
>    #   0.2974595   0.2638657
>
> I have spelled this out in detail to demonstrate that the key
> to accessing information in objects constructed by R lies in
> its structures (especially lists, vectors and matrices). You
> can find out what is involved for any function by looking for
> the section "Value" in its help page. For instance, the function
> summary() when applied to a GLM uses the method summary.glm(),
> so you can enter the command
>
>    ?summary.glm
>
> and then read what is in the section "Value". This shows that
> it is a list with components whose names are
>
>    call, family, deviance, ... , coefficients, ... , symbolic.cor
>
> and a component with name "Name" can be accessed using "$Name" as
> in "GLM$coef" (you can use "coef" instead of "coefficients" since
> the first four letters are [more than] enough to identify the name
> uniquely).
>
I would just add one suggestion to Ted's excellent tutorial:
R has the extractor function(s) coef() for getting the coefficients
(and SEs) for various types of models.

coef(GLM)
coef(summary(GLM))

While these will produce precisely the same output in the above
example, they may be the better way to go with, say, nonlinear
models. Using the first example in ?nls:

DNase1 <- subset(DNase, Run == 1)
fm1DNase1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1)

fm1DNase1$coef
# NULL  # <- probably not what was expected

coef(fm1DNase1)
#     Asym     xmid     scal
# 2.345180 1.483090 1.041455

Of course, looking at str(fm1DNase1) would show that there is no
component called "coefficients", but it might take a bit of head
scratching to realize that component "m" has as a subcomponent
the getAllPars() function which produces the ouput given by
coef(fm1DNase1).

I would recommend using extractor funtions like coef(), resid(),
etc. where available.

   -Peter Ehlers

> Once you get used to this, things become straightforward!
> Ted.
>
> --------------------------------------------------------------------
> E-Mail: (Ted Harding)<Ted.Harding at manchester.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 27-Feb-10                                       Time: 08:38:41
> ------------------------------ XFMail ------------------------------
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

-- 
Peter Ehlers
University of Calgary