[R] pglm package: fitted values and residuals

516 views
Skip to first unread message

alfonso...@uniparthenope.it

unread,
Apr 24, 2013, 4:11:15 AM4/24/13
to r-h...@r-project.org, yves.cr...@let.ish-lyon.cnrs.fr, Giovann...@generali.com
I'm using the package pglm and I'have estimated a "random probit model".
I need to save in a vector the fitted values and the residuals of the
model but I can not do it.

I tried with the command fitted.values using the following procedure
without results:

library(pglm)

m1_S<-pglm(Feed ~ Cons_PC_1 + imp_gen_1 + LGDP_PC_1 + lnEI_1 +
SH_Ren_1,data,family=binomial(probit),model="random",method="bfgs",index=c("Year","IDCountry"))

m1_S$fitted.values
residuals(m1)


Can someone help me about it?

Thanks


******************************************************************************
IL MERITO DEGLI STUDENTI VIENE RICONOSCIUTO

Il 5 per mille all'Universita' degli Studi di Napoli "Parthenope" incrementa le borse di studio agli studenti - codice fiscale 80018240632
http://www.uniparthenope.it/index.php/5xmille

http://www.uniparthenope.it/index.php/it/avvisi-sito-di-ateneo/2943-la-parthenope-premia-il-tuo-voto-di-diploma-ed-il-tuo-imegno-con-i-proventi-del-5-per-mille

Questa informativa e' inserita in automatico dal sistema al fine esclusivo della realizzazione dei fini istituzionali dell'ente.

______________________________________________
R-h...@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.

Paul Johnson

unread,
Apr 24, 2013, 3:39:21 PM4/24/13
to alfonso...@uniparthenope.it, R-help, yves.cr...@let.ish-lyon.cnrs.fr, Giovann...@generali.com
On Wed, Apr 24, 2013 at 3:11 AM, <alfonso...@uniparthenope.it> wrote:
> I'm using the package pglm and I'have estimated a "random probit model".
> I need to save in a vector the fitted values and the residuals of the model
> but I can not do it.
>
> I tried with the command fitted.values using the following procedure without
> results:
>
This is one of those "ask the pglm authors" questions. You should take
it up with the authors of the package. There is a specialized email
list R-sig-mixed where you will find more people working on this exact
same thing.

pglm looks like fun to me, but it is not quite done, so far as I can
tell. Well, the authors have not gone the "extra step" to make their
regression objects behave like other R regression objects. In case
you need alternative software, ask in R-sig-mixed. You'll learn that
most of these can be estimated with other packages. But I really like
the homogeneous user interface that is spelled out in pglm, and I
expect my students will run into the same questions that you have..

I just downloaded their source code, you probably ought to do that so
you can understand what they are doing. They provide the fitting
functions, but they do not do any of the other work necessary to make
these functions fit together with the R class framework. There are no
methods for "predict", anova, and so forth.

I'm in their R folder looking for implementations:

pauljohn@pols110:/tmp/pglm/R$ grep summary *
pauljohn@pols110:/tmp/pglm/R$ grep predict *
pauljohn@pols110:/tmp/pglm/R$ grep class *
pauljohn@pols110:/tmp/pglm/R$ grep fitted *
pglm.R: # glm models can be fitted

Run

> example(pglm)

what can we do after that?

> plot(anb)
Error in xy.coords(x, y, xlabel, ylabel, log) :
'x' is a list, but does not have components 'x' and 'y'

## Nothing.
## We do get a regression summary object, that's better than some
packages provide:

> anbsum <- summary(anb)

## And a coefficient table

> coef(anbsum)
Estimate Std. error t value Pr(> t)
(Intercept) -6.933764e-01 0.061391429 -11.294351205 1.399336e-29
wage 1.517009e-02 0.006375966 2.379261231 1.734738e-02
exper 1.314229e-03 0.007400129 0.177595444 8.590407e-01
ruralyes -8.594328e-05 0.051334716 -0.001674175 9.986642e-01

> model.matrix(anb)
Error in terms.default(object) : no terms component nor attribute
> anova(anb)
Error in UseMethod("anova") :
no applicable method for 'anova' applied to an object of class
"c('maxLik', 'maxim')"
> predict(anb)
Error in UseMethod("predict") :
no applicable method for 'predict' applied to an object of class
"c('maxLik', 'maxim')"

So, if you want those features with these models, you'll have to get
busy and do a lot of coding!

While working on regression support lately, I've reached the
conclusion that if an R package that claims to "do regression" but
does not provide methods for summary, predict, anova, nobs, fitted,
logLik, AIC, and so forth, then it is not done yet. Otherwise, users
like you who expect to be able to run methods like fitted or such have
a bad experience, as you are having now.

Maybe somebody reading this will remind us where the common list of R
regression methods is listed. I know for sure I've seen a document
about these things, but I'm baffled now trying to find it. But I'm
sure there is one.


pj

> library(pglm)
>
> m1_S<-pglm(Feed ~ Cons_PC_1 + imp_gen_1 + LGDP_PC_1 + lnEI_1 +
> SH_Ren_1,data,family=binomial(probit),model="random",method="bfgs",index=c("Year","IDCountry"))
>
> m1_S$fitted.values
> residuals(m1)
>
>
> Can someone help me about it?
>
> Thanks
>

Achim Zeileis

unread,
Apr 24, 2013, 5:37:07 PM4/24/13
to Paul Johnson, R-help, Giovann...@generali.com, yves.cr...@let.ish-lyon.cnrs.fr
On Wed, 24 Apr 2013, Paul Johnson wrote:

> On Wed, Apr 24, 2013 at 3:11 AM, <alfonso...@uniparthenope.it> wrote:
>
>> I'm using the package pglm and I'have estimated a "random probit
>> model". I need to save in a vector the fitted values and the residuals
>> of the model but I can not do it.
>>
>> I tried with the command fitted.values using the following procedure
>> without results:
>>
> This is one of those "ask the pglm authors" questions. You should take
> it up with the authors of the package. There is a specialized email
> list R-sig-mixed where you will find more people working on this exact
> same thing.
>
> pglm looks like fun to me, but it is not quite done, so far as I can
> tell. Well, the authors have not gone the "extra step" to make their
> regression objects behave like other R regression objects. In case you
> need alternative software, ask in R-sig-mixed. You'll learn that most of
> these can be estimated with other packages. But I really like the
> homogeneous user interface that is spelled out in pglm, and I expect my
> students will run into the same questions that you have..
>
> I just downloaded their source code, you probably ought to do that so
> you can understand what they are doing. They provide the fitting
> functions, but they do not do any of the other work necessary to make
> these functions fit together with the R class framework. There are no
> methods for "predict", anova, and so forth.

This is only partially true. In fact, "pglm" employs the framework
provided by the "maxLik" (by Ott Toomet and Arne Henningsen) and hence it
inherits some of the methods that "maxLik" provides for all of its fitted
model objects. So there is summary(), coef(), vcov(), AIC() work and you
can leverage tools like coeftest() from "lmtest", linearHypothesis() from
"car" or sandwich() from "sandwich" do work.

But it is certainly true that even more features would be desirable and I
think that Yves always planned on enhancing "pglm" at some point. (After
all it's still the initial version 0.1-0 on CRAN...)
I'm sure that there are many. One of my attempts to write up a list is in
Table 1 of vignette("betareg", package = "betareg").

Personally, I don't write anova() methods for my model objects because I
can leverage lrtest() and waldtest() from "lmtest" and linearHypothesis()
and deltaMethod() from "car" as long as certain standard methods are
available, including coef(), vcov(), logLik(), etc.

Similarly, an AIC() method is typically not needed as long as logLik() is
available. And BIC() works if nobs() is available in addition.

Best,
Z

Paul Johnson

unread,
Apr 25, 2013, 3:14:25 PM4/25/13
to Achim Zeileis, R-help, Giovann...@generali.com, yves.cr...@let.ish-lyon.cnrs.fr
On Wed, Apr 24, 2013 at 4:37 PM, Achim Zeileis <Achim....@uibk.ac.at> wrote:
> On Wed, 24 Apr 2013, Paul Johnson wrote:
>
>> On Wed, Apr 24, 2013 at 3:11 AM, <alfonso...@uniparthenope.it>
>> wrote:
>>
>>> I'm using the package pglm and I'have estimated a "random probit model".
>>> I need to save in a vector the fitted values and the residuals of the model
>>> but I can not do it.
>>>
>>> I tried with the command fitted.values using the following procedure
>>> without results:
>>>
>> This is one of those "ask the pglm authors" questions. You should take it
>> up with the authors of the package. There is a specialized email list
>> R-sig-mixed where you will find more people working on this exact same
>> thing.
>>
>> pglm looks like fun to me, but it is not quite done, so far as I can tell.
>
> I'm sure that there are many. One of my attempts to write up a list is in
> Table 1 of vignette("betareg", package = "betareg").

Yes! That's exactly the list I was thinking of. It was driving me
crazy I could not find it.

Thanks for the explanation. I don't think I should have implied that
the pglm author must actually implement all the methods, it is
certainly acceptable to leverage the methods that exist. It just
happened that the ones I tested were not implemented by any of the
affiliated packages.

But this thread leads me to one question I've wondered about recently.

Suppose I run somebody's regression function and out comes an object.

Do we have a way to ask that object "what are all of the methods that
might apply to you?" Here's why I wondered. You've noticed that
predict.lm has the interval="confidence" argument, but predict.glm
does not. So if I receive a regression model, I'd like to say to it
"do you have a predict method" and if I could get that predict method,
I could check to see if there is a formal argument interval. If it
does not, maybe I'd craft one for them.

pj



> Personally, I don't write anova() methods for my model objects because I can
> leverage lrtest() and waldtest() from "lmtest" and linearHypothesis() and
> deltaMethod() from "car" as long as certain standard methods are available,
> including coef(), vcov(), logLik(), etc.
>
> Similarly, an AIC() method is typically not needed as long as logLik() is
> available. And BIC() works if nobs() is available in addition.
>
> Best,
> Z
>
>>
>> pj
>>
>>> library(pglm)
>>>
>>> m1_S<-pglm(Feed ~ Cons_PC_1 + imp_gen_1 + LGDP_PC_1 + lnEI_1 +
>>>
>>> SH_Ren_1,data,family=binomial(probit),model="random",method="bfgs",index=c("Year","IDCountry"))
>>>
>>> m1_S$fitted.values
>>> residuals(m1)
>>>
>>>
>>> Can someone help me about it?
>>>
>>> Thanks
>>>
>>
>> ______________________________________________
>> R-h...@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.
>>
>



--
Paul E. Johnson
Professor, Political Science Assoc. Director
1541 Lilac Lane, Room 504 Center for Research Methods
University of Kansas University of Kansas
http://pj.freefaculty.org http://quant.ku.edu

Ista Zahn

unread,
Apr 25, 2013, 3:59:39 PM4/25/13
to Paul Johnson, R-help, Giovann...@generali.com, Achim Zeileis, yves.cr...@let.ish-lyon.cnrs.fr
Yes, minus the "might":

library(pglm)
example(pglm) # produces an object named "la"
sapply(class(la), function(x) methods(class=x)) # lists functions with
methods for objects of this class

Best,
Ista

Achim Zeileis

unread,
Apr 25, 2013, 4:45:36 PM4/25/13
to Ista Zahn, R-help, Giovann...@generali.com, yves.cr...@let.ish-lyon.cnrs.fr
Well, this shows you the methods that are available for the class but not
necessarily what arguments are supported. And even if the arguments are
available they do not necessarily mean the same thing. And some things may
or may not work via inheritance...

So coming back to Paul's question: Yes, I think it would be nice to have
support for this and in fact I have thought about similar infrastructure.
But so far I didn't have a good idea for a sufficiently robust/reliable
implementation. There are just so many details in the different model
objects that can be handled differently.

Best,
Z
Reply all
Reply to author
Forward
0 new messages