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.hard...@manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 27-Feb-10 Time: 08:38:41
------------------------------ XFMail ------------------------------
______________________________________________
R-help@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
______________________________________________
R-help@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.
