1# panelmodel and plm methods :
2
3## panelmodel methods :
4# - terms
5# - vcov
6# - fitted
7# - residuals
8# - df.residual
9# - coef
10# - print
11# - update
12# - deviance
13# - nobs
14
15## plm methods :
16# - summary
17# - print.summary
18# - predict
19# - formula
20# - plot
21# - residuals
22# - fitted
23
24
25#' @rdname plm
26#' @export
27terms.panelmodel <- function(x, ...){
28  terms(formula(x))
29}
30
31#' @rdname plm
32#' @export
33vcov.panelmodel <- function(object, ...){
34  object$vcov
35}
36
37#' @rdname plm
38#' @export
39fitted.panelmodel <- function(object, ...){
40  object$fitted.values
41}
42
43#' @rdname plm
44#' @export
45residuals.panelmodel <- function(object, ...){
46  object$residuals
47}
48
49#' @rdname plm
50#' @export
51df.residual.panelmodel <- function(object, ...){
52  object$df.residual
53}
54
55#' @rdname plm
56#' @export
57coef.panelmodel <- function(object, ...){
58  object$coefficients
59}
60
61#' @rdname plm
62#' @export
63print.panelmodel <- function(x, digits = max(3, getOption("digits") - 2),
64                             width = getOption("width"), ...){
65  cat("\nModel Formula: ")
66  print(formula(x))
67  cat("\nCoefficients:\n")
68  print(coef(x), digits = digits)
69  cat("\n")
70  invisible(x)
71}
72
73
74#' Extract Total Number of Observations Used in Estimated Panelmodel
75#'
76#' This function extracts the total number of 'observations' from a
77#' fitted panel model.
78#'
79#' The number of observations is usually the length of the residuals
80#' vector. Thus, `nobs` gives the number of observations actually
81#' used by the estimation procedure. It is not necessarily the number
82#' of observations of the model frame (number of rows in the model
83#' frame), because sometimes the model frame is further reduced by the
84#' estimation procedure. This is e.g. the case for first--difference
85#' models estimated by `plm(..., model = "fd")` where the model
86#' frame does not yet contain the differences (see also
87#' **Examples**).
88#'
89#' @name nobs.plm
90#' @aliases nobs
91#' @importFrom stats nobs
92#' @export nobs
93#' @param object a `panelmodel` object for which the number of
94#'     total observations is to be extracted,
95#' @param \dots further arguments.
96#' @return A single number, normally an integer.
97#' @seealso [pdim()]
98#' @keywords attribute
99#' @examples
100#'
101#' # estimate a panelmodel
102#' data("Produc", package = "plm")
103#' z <- plm(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp,data=Produc,
104#'          model="random", subset = gsp > 5000)
105#'
106#' nobs(z)       # total observations used in estimation
107#' pdim(z)$nT$N  # same information
108#' pdim(z)       # more information about the dimensions (no. of individuals and time periods)
109#'
110#' # illustrate difference between nobs and pdim for first-difference model
111#' data("Grunfeld", package = "plm")
112#' fdmod <- plm(inv ~ value + capital, data = Grunfeld, model = "fd")
113#' nobs(fdmod)      # 190
114#' pdim(fdmod)$nT$N # 200
115#'
116NULL
117
118# nobs() function to extract total number of observations used for estimating the panelmodel
119# like stats::nobs for lm objects
120# NB: here, use object$residuals rather than residuals(object)
121#     [b/c the latter could do NA padding once NA padding works for plm objects.
122#      NA padded residuals would yield wrong result for nobs!]
123
124#' @rdname nobs.plm
125#' @export
126nobs.panelmodel <- function(object, ...) {
127  if (inherits(object, "plm") || inherits(object, "panelmodel")) return(length(object$residuals))
128    else stop("Input 'object' needs to be of class 'plm' or 'panelmodel'")
129}
130
131# No of obs calculated as in print.summary.pgmm [code copied from there]
132#' @rdname nobs.plm
133#' @export
134nobs.pgmm <- function(object, ...) {
135  if (inherits(object, "pgmm")) return(sum(unlist(object$residuals, use.names = FALSE) != 0))
136    else stop("Input 'object' needs to be of class 'pgmm', i. e., a GMM estimation with panel data estimated by pgmm()")
137}
138
139
140
141
142# Almost the same as the default method except that update.formula is
143# replaced by update, so that the Formula method is used to update the
144# formula
145
146#' @rdname plm
147#' @export
148update.panelmodel <- function (object, formula., ..., evaluate = TRUE){
149    if (is.null(call <- object$call)) # was: getCall(object)))
150        stop("need an object with call component")
151    extras <- match.call(expand.dots = FALSE)$...
152    # update.Formula fails if latter rhs are . ; simplify the formula
153    # by removing the latter parts
154
155    if (! missing(formula.)){
156        newform <- Formula(formula.)
157        if (length(newform)[2L] == 2L && attr(newform, "rhs")[2L] == as.name("."))
158            newform <- formula(newform, rhs = 1)
159        call$formula <- update(formula(object), newform)
160    }
161    if (length(extras)) {
162        existing <- !is.na(match(names(extras), names(call)))
163        for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
164        if (any(!existing)) {
165            call <- c(as.list(call), extras[!existing])
166            call <- as.call(call)
167        }
168    }
169    if (evaluate)
170        eval(call, parent.frame())
171    else call
172}
173
174#' @rdname plm
175#' @export
176deviance.panelmodel <- function(object, model = NULL, ...){
177  if (is.null(model)) as.numeric(crossprod(resid(object)))
178  else as.numeric(crossprod(residuals(object, model = model)))
179}
180
181
182
183# summary.plm creates a specific summary.plm object that is derived
184# from the associated plm object
185
186
187#' Summary for plm objects
188#'
189#' The summary method for plm objects generates some more information about
190#' estimated plm models.
191#'
192#' The `summary` method for plm objects (`summary.plm`) creates an
193#' object of class `c("summary.plm", "plm", "panelmodel")` that
194#' extends the plm object it is run on with various information about
195#' the estimated model like (inferential) statistics, see
196#' **Value**. It has an associated print method
197#' (`print.summary.plm`).
198#'
199#' @aliases summary.plm
200#' @param object an object of class `"plm"`,
201#' @param x an object of class `"summary.plm"`,
202#' @param subset a character or numeric vector indicating a subset of
203#'     the table of coefficients to be printed for
204#'     `"print.summary.plm"`,
205#' @param vcov a variance--covariance matrix furnished by the user or
206#'     a function to calculate one (see **Examples**),
207#' @param digits number of digits for printed output,
208#' @param width the maximum length of the lines in the printed output,
209#' @param eq the selected equation for list objects
210#' @param \dots further arguments.
211#' @return An object of class `c("summary.plm", "plm",
212#'     "panelmodel")`.  Some of its elements are carried over from the
213#'     associated plm object and described there
214#'     ([plm()]). The following elements are new or changed
215#'     relative to the elements of a plm object:
216#'
217#' \item{fstatistic}{'htest' object: joint test of significance of
218#' coefficients (F or Chi-square test) (robust statistic in case of
219#' supplied argument `vcov`, see [pwaldtest()] for details),}
220#'
221#' \item{coefficients}{a matrix with the estimated coefficients,
222#' standard errors, t--values, and p--values, if argument `vcov` was
223#' set to non-`NULL` the standard errors (and t-- and p--values) in
224#' their respective robust variant,}
225#'
226#' \item{vcov}{the "regular" variance--covariance matrix of the coefficients (class "matrix"),}
227#'
228#' \item{rvcov}{only present if argument `vcov` was set to non-`NULL`:
229#' the furnished variance--covariance matrix of the coefficients
230#' (class "matrix"),}
231#'
232#' \item{r.squared}{a named numeric containing the R-squared ("rsq")
233#' and the adjusted R-squared ("adjrsq") of the model,}
234#'
235#' \item{df}{an integer vector with 3 components, (p, n-p, p*), where
236#' p is the number of estimated (non-aliased) coefficients of the
237#' model, n-p are the residual degrees of freedom (n being number of
238#' observations), and p* is the total number of coefficients
239#' (incl. any aliased ones).}
240#'
241#' @export
242#' @author Yves Croissant
243#' @seealso [plm()] for estimation of various models; [vcovHC()] for
244#'     an example of a robust estimation of variance--covariance
245#'     matrix; [r.squared()] for the function to calculate R-squared;
246#'     [stats::print.power.htest()] for some information about class
247#'     "htest"; [fixef()] to compute the fixed effects for "within"
248#'     (=fixed effects) models and [within_intercept()] for an
249#'     "overall intercept" for such models; [pwaldtest()]
250#' @keywords regression
251#' @examples
252#'
253#' data("Produc", package = "plm")
254#' zz <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,
255#'           data = Produc, index = c("state","year"))
256#' summary(zz)
257#'
258#' # summary with a furnished vcov, passed as matrix, as function, and
259#' # as function with additional argument
260#' data("Grunfeld", package = "plm")
261#' wi <- plm(inv ~ value + capital,
262#'           data = Grunfeld, model="within", effect = "individual")
263#' summary(wi, vcov = vcovHC(wi))
264#' summary(wi, vcov = vcovHC)
265#' summary(wi, vcov = function(x) vcovHC(x, method = "white2"))
266#'
267#' # extract F statistic
268#' wi_summary <- summary(wi)
269#' Fstat <- wi_summary[["fstatistic"]]
270#'
271#' # extract estimates and p-values
272#' est <- wi_summary[["coefficients"]][ , "Estimate"]
273#' pval <- wi_summary[["coefficients"]][ , "Pr(>|t|)"]
274#'
275#' # print summary only for coefficent "value"
276#' print(wi_summary, subset = "value")
277#'
278summary.plm <- function(object, vcov = NULL, ...){
279
280  vcov_arg <- vcov
281  model <- describe(object, "model")
282  effect <- describe(object, "effect")
283  random.method <- describe(object, "random.method")
284
285  # determine if intercept-only model (no other regressors)
286  coef_wo_int <- object$coefficients[!(names(coef(object)) %in% "(Intercept)")]
287  int.only <- !length(coef_wo_int)
288
289  # as cor() is not defined for intercept-only models, use different approach
290  # for R-squared ("rss" and "ess" are defined)
291  object$r.squared <- if(!int.only) {
292      c(rsq    = r.squared(object),
293        adjrsq = r.squared(object, dfcor = TRUE))
294    } else {
295      c(rsq    = r.squared(object, type = "rss"),
296        adjrsq = r.squared(object, type = "rss", dfcor = TRUE))
297    }
298
299  ## determine if standard normal and Chisq test or t distribution and F test to be used
300  ## (normal/chisq for all random models, all IV models, and HT via plm(., model="ht"))
301  use.norm.chisq <- if(model == "random" ||
302                       length(formula(object))[2L] >= 2L ||
303                       model == "ht") TRUE else FALSE
304
305  # perform Wald test of joint sign. of regressors only if there are
306  # other regressors besides the intercept
307  if(!int.only) {
308    object$fstatistic <- pwaldtest(object,
309                                   test = if(use.norm.chisq) "Chisq" else "F",
310                                   vcov = vcov_arg)
311  }
312
313
314  # construct the table of coefficients
315  if (!is.null(vcov_arg)) {
316    if (is.matrix(vcov_arg))   rvcov <- vcov_arg
317    if (is.function(vcov_arg)) rvcov <- vcov_arg(object)
318    std.err <- sqrt(diag(rvcov))
319  } else {
320    std.err <- sqrt(diag(stats::vcov(object)))
321  }
322  b <- coefficients(object)
323  z <- b / std.err
324  p <- if(use.norm.chisq) {
325    2 * pnorm(abs(z), lower.tail = FALSE)
326  } else {
327    2 * pt(abs(z), df = object$df.residual, lower.tail = FALSE)
328  }
329
330  # construct the object of class summary.plm
331  object$coefficients <- cbind(b, std.err, z, p)
332  colnames(object$coefficients) <- if(use.norm.chisq) {
333    c("Estimate", "Std. Error", "z-value", "Pr(>|z|)")
334  } else { c("Estimate", "Std. Error", "t-value", "Pr(>|t|)") }
335
336  ## add some info to summary.plm object
337  # robust vcov (next to "normal" vcov)
338  if (!is.null(vcov_arg)) {
339    object$rvcov <- rvcov
340    rvcov.name <- paste0(deparse(substitute(vcov)))
341    attr(object$rvcov, which = "rvcov.name") <- rvcov.name
342  }
343
344  # mimics summary.lm's 'df' component
345  # 1st entry: no. coefs (w/o aliased coefs); 2nd: residual df; 3rd no. coefs /w aliased coefs
346  # NB: do not use length(object$coefficients) for 3rd entry!
347  object$df <- c(length(b), object$df.residual, length(object$aliased))
348
349  class(object) <- c("summary.plm", "plm", "panelmodel")
350  object
351}
352
353#' @rdname summary.plm
354#' @export
355print.summary.plm <- function(x, digits = max(3, getOption("digits") - 2),
356                              width = getOption("width"), subset = NULL, ...){
357  formula <- formula(x)
358  has.instruments <- (length(formula)[2L] >= 2L)
359  effect <- describe(x, "effect")
360  model  <- describe(x, "model")
361  if (model != "pooling") { cat(paste(effect.plm.list[effect], " ", sep = "")) }
362  cat(paste(model.plm.list[model], " Model", sep = ""))
363
364  if (model == "random"){
365    ercomp <- describe(x, "random.method")
366    cat(paste(" \n   (",
367              random.method.list[ercomp],
368              "'s transformation)\n",
369              sep = ""))
370  }
371  else{
372    cat("\n")
373  }
374
375  if (has.instruments){
376    cat("Instrumental variable estimation\n")
377    if(model != "within") {
378      # don't print transformation method for FE models as there is only one
379      # such method for FE models but plenty for other model types
380      ivar <- describe(x, "inst.method")
381      cat(paste0("   (", inst.method.list[ivar], "'s transformation)\n"))
382    }
383  }
384
385  if (!is.null(x$rvcov)) {
386    cat("\nNote: Coefficient variance-covariance matrix supplied: ", attr(x$rvcov, which = "rvcov.name"), "\n", sep = "")
387  }
388
389  cat("\nCall:\n")
390  print(x$call)
391  cat("\n")
392  pdim <- pdim(x)
393  print(pdim)
394  if (model %in% c("fd", "between")) {
395    # print this extra info, b/c model.frames of FD and between models
396    # have original (undifferenced/"un-between-ed") obs/rows of the data
397    cat(paste0("Observations used in estimation: ", nobs(x), "\n"))}
398
399  if (model == "random"){
400    cat("\nEffects:\n")
401    print(x$ercomp)
402  }
403  cat("\nResiduals:\n")
404  df <- x$df
405  rdf <- df[2L]
406  if (rdf > 5L) {
407    save.digits <- unlist(options(digits = digits))
408    on.exit(options(digits = save.digits))
409    print(sumres(x))
410  } else if (rdf > 0L) print(residuals(x), digits = digits)
411  if (rdf == 0L) { # estimation is a perfect fit
412    cat("ALL", x$df[1L], "residuals are 0: no residual degrees of freedom!")
413    cat("\n")
414  }
415
416  if (any(x$aliased, na.rm = TRUE)) {
417    # na.rm = TRUE because currently, RE tw unbalanced models might have NAs?
418    naliased <- sum(x$aliased, na.rm = TRUE)
419    cat("\nCoefficients: (", naliased, " dropped because of singularities)\n", sep = "")
420  } else cat("\nCoefficients:\n")
421
422  if (is.null(subset)) printCoefmat(coef(x), digits = digits)
423  else printCoefmat(coef(x)[subset, , drop = FALSE], digits = digits)
424  cat("\n")
425  cat(paste("Total Sum of Squares:    ", signif(tss(x),      digits), "\n", sep = ""))
426  cat(paste("Residual Sum of Squares: ", signif(deviance(x), digits), "\n", sep = ""))
427  cat(paste("R-Squared:      ", signif(x$r.squared[1L], digits),      "\n", sep = ""))
428  cat(paste("Adj. R-Squared: ", signif(x$r.squared[2L], digits),      "\n", sep = ""))
429
430  # print Wald test of joint sign. of regressors only if there is a statistic
431  # in summary.plm object (not computed by summary.plm if there are no other
432  # regressors than the intercept
433  if(!is.null(fstat <- x$fstatistic)) {
434    if (names(fstat$statistic) == "F"){
435      cat(paste("F-statistic: ", signif(fstat$statistic),
436                " on ", fstat$parameter["df1"]," and ", fstat$parameter["df2"],
437                " DF, p-value: ", format.pval(fstat$p.value,digits=digits), "\n", sep=""))
438    }
439    else{
440      cat(paste("Chisq: ", signif(fstat$statistic),
441                " on ", fstat$parameter,
442                " DF, p-value: ", format.pval(fstat$p.value, digits = digits), "\n", sep=""))
443    }
444  }
445  invisible(x)
446}
447
448#' @rdname plm
449#' @export
450predict.plm <- function(object, newdata = NULL, ...){
451  tt <- terms(object)
452  if (is.null(newdata)){
453    result <- fitted(object, ...)
454  }
455  else{
456    Terms <- delete.response(tt)
457    m <- model.frame(Terms, newdata)
458    X <- model.matrix(Terms, m)
459    beta <- coef(object)
460    result <- as.numeric(crossprod(beta, t(X)))
461  }
462  result
463}
464
465#' @rdname plm
466#' @export
467formula.plm <- function(x, ...){
468  x$formula
469}
470
471#' @rdname plm
472#' @export
473plot.plm <- function(x, dx = 0.2, N = NULL, seed = 1,
474                     within = TRUE, pooling = TRUE,
475                     between = FALSE, random = FALSE, ...){
476    set.seed(seed)# 8 est bien pour beertax
477    subs <- ! is.null(N)
478    x <- update(x, model = "within")
479    mco <- update(x, model = "pooling")
480    if (random) re <- update(x, model = "random")
481    if (between) be <- update(x, model = "between")
482    pdim <- pdim(x)
483    n <- pdim$nT$n
484    if (! subs) N <- n
485    ids <- unique(index(x, "id"))
486    if (subs) ids <- ids[sample(1:length(ids), N, replace = FALSE)]
487    sel <- index(x, "id") %in% ids
488    T. <- pdim$nT$T
489    cols <- rainbow(N)
490    pts <- sample(1:25, N, replace = TRUE)
491    thex <- as.numeric(model.matrix(x, model = "pooling")[sel, 2L])
492    they <- as.numeric(pmodel.response(x, model = "pooling")[sel])
493    plot(thex, they, col = rep(cols, each = T.),
494         pch = rep(pts, each = T.), ann = FALSE, las = 1)
495    idsel <- as.numeric(index(x, "id")[sel])
496    meanx <- tapply(thex, idsel, mean)
497    meany <- tapply(they, idsel, mean)
498    points(meanx, meany, pch = 19, col = cols, cex = 1.5)
499    if (within){
500        beta <- coef(x)
501        alphas <- meany - meanx * beta
502        dx <- dx * (max(thex) - min(thex))
503        for (i in 1:N){
504            xmin <- meanx[i] - dx
505            xmax <- meanx[i] + dx
506            ymin <- alphas[i] + beta * xmin
507            ymax <- alphas[i] + beta * xmax
508            lines(c(xmin, xmax), c(ymin, ymax), col = cols[i])
509        }
510    }
511    if(random) abline(coef(re)[1L], coef(re)[2L], lty = "dotted")
512    if(pooling) abline(coef(mco), lty = "dashed")
513    if(between) abline(coef(be), lty = "dotdash")
514    # where to put the legends, depends on the sign of the OLS slope
515    modploted <- c(random, pooling, between, within)
516    if (sum(modploted)){
517        poslegend <- ifelse(beta > 0, "topleft", "topright")
518        ltylegend <- c("dotted", "dashed", "dotdash", "solid")[modploted]
519        leglegend <- c("random", "pooling", "between", "within")[modploted]
520        legend(poslegend, lty = ltylegend, legend = leglegend)
521    }
522}
523
524#' @rdname plm
525#' @export
526residuals.plm <- function(object, model = NULL, effect = NULL,  ...){
527    if (is.null(model) && is.null(effect)){
528        model <- describe(object, "model")
529        res <- object$residuals
530    }
531    else{
532        cl <- match.call(expand.dots = FALSE)
533        # fitted -> call to the plm method, used to be fitted.plm
534        # which is not exported
535#        cl[[1L]] <- as.name("fitted.plm")
536        cl[[1L]] <- as.name("fitted")
537        bX <- eval(cl, parent.frame())
538        if (is.null(model))  model  <- describe(object, "model")
539        if (is.null(effect)) effect <- describe(object, "effect")
540        y <- pmodel.response(object, model = model, effect = effect)
541        res <- y - bX
542    }
543    res <- if (model %in% c("between", "fd")) {
544      # these models "compress" the data, thus an index does not make sense here
545      # -> do not return pseries but plain numeric
546      res
547    } else {
548      structure(res, index = index(object), class = union("pseries", class(res)))
549    }
550    return(res)
551}
552
553#' @rdname plm
554#' @export
555fitted.plm <- function(object, model = NULL, effect = NULL, ...){
556    fittedmodel <- describe(object, "model")
557    if (is.null(model)) model <- fittedmodel
558    if (is.null(effect)) effect <- describe(object, "effect")
559    if (fittedmodel == "random") theta <- ercomp(object)$theta else theta <- NULL
560    X <- model.matrix(object, model = "pooling")
561    y <- pmodel.response(object, model = "pooling", effect = effect)
562    beta <- coef(object)
563    comonpars <- intersect(names(beta), colnames(X))
564    bX <- as.numeric(crossprod(t(X[, comonpars, drop = FALSE]), beta[comonpars]))
565    bX <- structure(bX, index = index(object), class = union("pseries", class(bX)))
566    if (fittedmodel == "within"){
567        intercept <- mean(y - bX)
568        bX <- bX + intercept
569    }
570    ptransform(bX, model = model, effect = effect, theta = theta)
571}
572