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