1 2#' Hausman Test for Panel Models 3#' 4#' Specification test for panel models. 5#' 6#' The Hausman test (sometimes also called Durbin--Wu--Hausman test) 7#' is based on the difference of the vectors of coefficients of two 8#' different models. The `panelmodel` method computes the original 9#' version of the test based on a quadratic form 10#' \insertCite{HAUS:78}{plm}. The `formula` method, if 11#' `method="chisq"` (default), computes the original version of the 12#' test based on a quadratic form; if `method="aux"` then the 13#' auxiliary-regression-based version in Wooldridge (2010, 14#' Sec. 10.7.3.) is computed instead \insertCite{@WOOL:10 Sec.10.7.3}{plm}. 15#' Only the latter can be robustified by specifying a robust 16#' covariance estimator as a function through the argument `vcov` (see 17#' **Examples**). 18#' 19#' The equivalent tests in the **one-way** case using a between 20#' model (either "within vs. between" or "random vs. between") 21#' \insertCite{@see @HAUS:TAYL:81 or @BALT:13 Sec.4.3}{plm} can also 22#' be performed by `phtest`, but only for `test = "chisq"`, not for 23#' the regression-based test. NB: These equivalent tests using the 24#' between model do not extend to the two-ways case. There are, 25#' however, some other equivalent tests, 26#' \insertCite{@see @KANG:85 or @BALT:13 Sec.4.3.7}{plm}, 27#' but those are unsupported by `phtest`. 28#' 29#' @aliases phtest 30#' @param x an object of class `"panelmodel"` or `"formula"`, 31#' @param x2 an object of class `"panelmodel"`, 32#' @param model a character vector containing the names of two models 33#' (length(model) must be 2), 34#' @param data a `data.frame`, 35#' @param method one of `"chisq"` or `"aux"`, 36#' @param index an optional vector of index variables, 37#' @param vcov an optional covariance function, 38#' @param \dots further arguments to be passed on. For the formula method, 39#' place argument `effect` here to compare, e.g., twoway models 40#' (`effect = "twoways"`) Note: Argument `effect` is not respected in 41#' the panelmodel method. 42#' @return An object of class `"htest"`. 43#' @export 44#' @author Yves Croissant, Giovanni Millo 45#' @references 46#' 47#' \insertRef{HAUS:78}{plm} 48#' 49#' \insertRef{HAUS:TAYL:81}{plm} 50#' 51#' \insertRef{KANG:85}{plm} 52#' 53#' \insertRef{WOOL:10}{plm} 54#' 55#' \insertRef{BALT:13}{plm} 56#' 57#' @keywords htest 58#' @examples 59#' 60#' data("Gasoline", package = "plm") 61#' form <- lgaspcar ~ lincomep + lrpmg + lcarpcap 62#' wi <- plm(form, data = Gasoline, model = "within") 63#' re <- plm(form, data = Gasoline, model = "random") 64#' phtest(wi, re) 65#' phtest(form, data = Gasoline) 66#' phtest(form, data = Gasoline, method = "aux") 67#' 68#' # robust Hausman test (regression-based) 69#' phtest(form, data = Gasoline, method = "aux", vcov = vcovHC) 70#' 71#' # robust Hausman test with vcov supplied as a 72#' # function and additional parameters 73#' phtest(form, data = Gasoline, method = "aux", 74#' vcov = function(x) vcovHC(x, method="white2", type="HC3")) 75#' 76phtest <- function(x,...){ 77 UseMethod("phtest") 78} 79 80#' @rdname phtest 81#' @export 82phtest.formula <- function(x, data, model = c("within", "random"), 83 method = c("chisq", "aux"), 84 index = NULL, vcov = NULL, ...) { 85 # TODO: No argument 'effect' here, maybe introduce? 86 # it gets evaluated due to the eval() call for method="chisq" 87 # and since rev. 305 due to extraction from dots (...) in method="aux" as a quick fix 88 # If introduced as argument, change doc accordingly (currently, effect arg is mentioned in ...) 89 90 if (length(model) != 2) stop("two models should be indicated in argument 'model'") 91 for (i in 1:2){ 92 model.name <- model[i] 93 if(!(model.name %in% names(model.plm.list))){ 94 stop("model must be one of ", oneof(model.plm.list)) 95 } 96 } 97 switch(match.arg(method), 98 chisq={ 99 cl <- match.call(expand.dots = TRUE) 100 cl$model <- model[1L] 101 names(cl)[2L] <- "formula" 102 m <- match(plm.arg, names(cl), 0L) 103 cl <- cl[c(1L, m)] 104 cl[[1L]] <- as.name("plm") 105 plm.model.1 <- eval(cl, parent.frame()) 106 plm.model.2 <- update(plm.model.1, model = model[2L]) 107 return(phtest(plm.model.1, plm.model.2)) 108 }, 109 aux={ 110 ## some interface checks here 111 if (model[1L] != "within") { 112 stop("Please supply 'within' as first model type") 113 } 114 115 if (!is.null(vcov) && !is.function(vcov)) stop("argument 'vcov' needs to be a function") 116 117 ## set pdata.frame 118 if (!inherits(data, "pdata.frame")) data <- pdata.frame(data, index = index) #, ...) 119 120 row.names(data) <- NULL # reset rownames of original data set (->numbers rownames in clean sequence) to make rownames 121 # comparable for later comparison to obs used in estimation of models (get rid of NA values) 122 # [needed because pmodel.response() and model.matrix() do not retain fancy rownames, but rownames] 123 124 # rev. 305: quick and dirty fix for missing effect argument in function 125 # signature for formula interface/test="aux": see if effect is in dots and extract 126 dots <- list(...) 127 # print(dots) # DEBUG printing 128 effect <- if(!is.null(dots$effect)) dots$effect else NULL 129 # calculate FE and RE model 130 131 fe_mod <- plm(formula = x, data = data, model = model[1L], effect = effect) 132 re_mod <- plm(formula = x, data = data, model = model[2L], effect = effect) 133 134 ## DEBUG printing: 135 # print(effect) 136 # print(model) 137 # print(paste0("mod1: ", describe(fe_mod, "effect"))) 138 # print(paste0("mod2: ", describe(re_mod, "effect"))) 139 # print(fe_mod) 140 # print(re_mod) 141 reY <- pmodel.response(re_mod) 142# reX <- model.matrix(re_mod)[ , -1, drop = FALSE] # intercept not needed; drop=F needed to prevent matrix 143# feX <- model.matrix(fe_mod, cstcovar.rm = TRUE) # from degenerating to vector if only one regressor 144 reX <- model.matrix(re_mod, cstcovar.rm = "intercept") 145 feX <- model.matrix(fe_mod, cstcovar.rm = "all") 146 147 dimnames(feX)[[2L]] <- paste(dimnames(feX)[[2L]], "tilde", sep=".") 148 ## estimated models could have fewer obs (due dropping of NAs) compared to the original data 149 ## => match original data and observations used in estimated models 150 ## routine adapted from lmtest::bptest 151 commonrownames <- intersect(intersect(intersect(row.names(data), names(reY)), row.names(reX)), row.names(feX)) 152 if (!(all(c(row.names(data) %in% commonrownames, commonrownames %in% row.names(data))))) { 153 data <- data[commonrownames, ] 154 reY <- reY[commonrownames] 155 reX <- reX[commonrownames, ] 156 feX <- feX[commonrownames, ] 157 } 158 159 # Tests of correct matching of obs (just for safety ...) 160 if (!all.equal(length(reY), nrow(data), nrow(reX), nrow(feX))) 161 stop("number of cases/observations do not match, most likely due to NAs in \"data\"") 162 if (any(c(is.na(names(reY)), is.na(row.names(data)), is.na(row.names(reX)), is.na(row.names(feX))))) 163 stop("one (or more) rowname(s) is (are) NA") 164 if (!all.equal(names(reY), row.names(data), row.names(reX), row.names(feX))) 165 stop("row.names of cases/observations do not match, most likely due to NAs in \"data\"") 166 167 ## fetch indices here, check pdata 168 ## construct data set and formula for auxiliary regression 169 data <- pdata.frame(cbind(index(data), reY, reX, feX)) 170 auxfm <- as.formula(paste("reY~", 171 paste(dimnames(reX)[[2L]], 172 collapse="+"), "+", 173 paste(dimnames(feX)[[2L]], 174 collapse="+"), sep="")) 175 auxmod <- plm(formula = auxfm, data = data, model = "pooling") 176 nvars <- dim(feX)[[2L]] 177 R <- diag(1, nvars) 178 r <- rep(0, nvars) # here just for clarity of illustration 179 omega0 <- vcov(auxmod)[(nvars+2):(nvars*2+1), 180 (nvars+2):(nvars*2+1)] 181 Rbr <- R %*% coef(auxmod)[(nvars+2):(nvars*2+1)] - r 182 183 h2t <- as.numeric(crossprod(Rbr, solve(omega0, Rbr))) 184 ph2t <- pchisq(h2t, df = nvars, lower.tail = FALSE) 185 186 df <- nvars 187 names(df) <- "df" 188 names(h2t) <- "chisq" 189 190 if (!is.null(vcov)) { 191 vcov <- paste(", vcov: ", 192 paste(deparse(substitute(vcov))), 193 sep="") 194 } 195 196 haus2 <- list(statistic = h2t, 197 p.value = ph2t, 198 parameter = df, 199 method = paste("Regression-based Hausman test", 200 vcov, sep=""), 201 alternative = "one model is inconsistent", 202 data.name = paste(deparse(substitute(x)))) 203 class(haus2) <- "htest" 204 return(haus2) 205 }) 206} 207 208#' @rdname phtest 209#' @export 210phtest.panelmodel <- function(x, x2, ...){ 211 coef.wi <- coef(x) 212 coef.re <- coef(x2) 213 vcov.wi <- vcov(x) 214 vcov.re <- vcov(x2) 215 names.wi <- names(coef.wi) 216 names.re <- names(coef.re) 217 common_coef_names <- names.re[names.re %in% names.wi] 218 coef.h <- common_coef_names[!(common_coef_names %in% "(Intercept)")] # drop intercept if included (relevant when between model input) 219 if(length(coef.h) == 0L) stop("no common coefficients in models") 220 dbeta <- coef.wi[coef.h] - coef.re[coef.h] 221 df <- length(dbeta) 222 dvcov <- vcov.wi[coef.h, coef.h] - vcov.re[coef.h, coef.h] 223 224 #### BEGIN cater for equivalent test within vs. between 225 # Baltagi (2013), Sec. 4.3, pp. 77, 81 226 modx <- describe(x, what = "model") 227 modx2 <- describe(x2, what = "model") 228 effx <- describe(x, what = "effect") 229 effx2 <- describe(x2, what = "effect") 230 231 # Tests with between model do not extend to two-ways case -> give error 232 # There are, however, some equiv. tests with the individual/time between 233 # model, but let's not support them (see Kang (1985), Baltagi (2013), Sec. 4.3.7) 234 if ( (modx == "between" || modx2 == "between") 235 && (effx == "twoways" || effx2 == "twoways")) stop("tests with between model in twoways case not supported") 236 237 # in case of one-way within vs. between (m3 in Baltagi (2013), pp. 77, 81) 238 # the variances need to be added (not subtracted like in the other cases) 239 if ( (modx == "within" && modx2 == "between") 240 || (modx2 == "within" && modx == "between")) { 241 dvcov <- vcov.wi[coef.h, coef.h] + vcov.re[coef.h, coef.h] 242 } 243 #### END cater for equivalent tests with between model 244 245 stat <- as.numeric(abs(t(dbeta) %*% solve(dvcov) %*% dbeta)) 246 pval <- pchisq(stat, df = df, lower.tail = FALSE) 247 names(stat) <- "chisq" 248 parameter <- df 249 names(parameter) <- "df" 250 alternative <- "one model is inconsistent" 251 252 ## DEBUG printing: 253 # print(paste0("mod1: ", describe(x, "effect"))) 254 # print(paste0("mod2: ", describe(x2, "effect"))) 255 256 res <- list(statistic = stat, 257 p.value = pval, 258 parameter = parameter, 259 method = "Hausman Test", 260 data.name = data.name(x), 261 alternative = alternative) 262 class(res) <- "htest" 263 return(res) 264} 265 266############## plmtest() ############################################ 267# For a concise overview with original references, see 268# Baltagi (2013), Econometric Analysis of Panel Data, 5th edition, pp. 68-76 (balanced), pp. 200-203 (unbalanced). 269# 270# balanced (original) version of Breusch-Pagan test: 271# T.S. Breusch & A.R. Pagan (1979), 272# A Simple Test for Heteroscedasticity and Random Coefficient Variation, 273# Econometrica 47, pp. 1287-1294 274# 275# unbalanced version: 276# Baltagi/Li (1990), 277# A lagrange multiplier test for the error components model with incomplete panels, 278# Econometric Reviews, 9, pp. 103-107, 279 280 281# pchibarsq: helper function: "p-function" for mixed chisq (also called chi-bar-squared) 282# used in plmtest(., type = "ghm"), see Baltagi (2013), pp. 71-72, 74, 88, 202-203, 209 283# 284# a reference for the distribution seems to be 285# Dykstra, R./El Barmi, H., Chi-Bar-Square Distributions, in: Encyclopedia of Statistical Sciences, 286# DOI: 10.1002/0471667196.ess0265.pub2 287pchibarsq <- function(q, df, weights, lower.tail = TRUE, ... ) { 288 # NB: other parameters in dots (...): not checked if valid! (ncp, log, ...) 289 res <- sum(weights * pchisq(q, df = df, lower.tail = lower.tail, ...)) 290 return(res) 291} 292 293 294 295 296#' Lagrange FF Multiplier Tests for Panel Models 297#' 298#' Test of individual and/or time effects for panel models. 299#' 300#' These Lagrange multiplier tests use only the residuals of the 301#' pooling model. The first argument of this function may be either a 302#' pooling model of class `plm` or an object of class `formula` 303#' describing the model. For inputted within (fixed effects) or random 304#' effects models, the corresponding pooling model is calculated 305#' internally first as the tests are based on the residuals of the 306#' pooling model. 307#' 308#' The `"bp"` test for unbalanced panels was derived in 309#' \insertCite{BALT:LI:90;textual}{plm} 310#' (1990), the `"kw"` test for unbalanced panels in 311#' \insertCite{BALT:CHAN:LI:98;textual}{plm}. 312#' 313#' The `"ghm"` test and the `"kw"` test were extended to two-way 314#' effects in \insertCite{BALT:CHAN:LI:92;textual}{plm}. 315#' 316#' For a concise overview of all these statistics see 317#' \insertCite{BALT:03;textual}{plm}, Sec. 4.2, pp. 68--76 (for balanced 318#' panels) and Sec. 9.5, pp. 200--203 (for unbalanced panels). 319#' 320#' @aliases plmtest 321#' @param x an object of class `"plm"` or a formula of class 322#' `"formula"`, 323#' @param data a `data.frame`, 324#' @param effect a character string indicating which effects are 325#' tested: individual effects (`"individual"`), time effects 326#' (`"time"`) or both (`"twoways"`), 327#' @param type a character string indicating the test to be performed: 328#' 329#' - `"honda"` (default) for \insertCite{HOND:85;textual}{plm}, 330#' - `"bp"` for \insertCite{BREU:PAGA:80;textual}{plm}, 331#' - `"kw"` for \insertCite{KING:WU:97;textual}{plm}, or 332#' - `"ghm"` for \insertCite{GOUR:HOLL:MONF:82;textual}{plm} for 333#' unbalanced panel data sets, the respective unbalanced version 334#' of the tests are computed, 335#' 336#' @param \dots further arguments passed to `plmtest`. 337#' @return An object of class `"htest"`. 338#' @note For the King-Wu statistics (`"kw"`), the oneway statistics 339#' (`"individual"` and `"time"`) coincide with the respective 340#' Honda statistics (`"honda"`); twoway statistics of `"kw"` and 341#' `"honda"` differ. 342#' @export 343#' @author Yves Croissant (initial implementation), Kevin Tappe 344#' (generalization to unbalanced panels) 345#' @seealso [pFtest()] for individual and/or time effects tests based 346#' on the within model. 347#' @references 348#' 349#' \insertRef{BALT:13}{plm} 350#' 351#' \insertRef{BALT:LI:90}{plm} 352#' 353#' \insertRef{BALT:CHAN:LI:92}{plm} 354#' 355#' \insertRef{BALT:CHAN:LI:98}{plm} 356#' 357#' \insertRef{BREU:PAGA:80}{plm} 358#' 359#' \insertRef{GOUR:HOLL:MONF:82}{plm} 360#' 361#' \insertRef{HOND:85}{plm} 362#' 363#' \insertRef{KING:WU:97}{plm} 364#' 365#' @keywords htest 366#' @examples 367#' 368#' data("Grunfeld", package = "plm") 369#' g <- plm(inv ~ value + capital, data = Grunfeld, model = "pooling") 370#' plmtest(g) 371#' plmtest(g, effect="time") 372#' plmtest(inv ~ value + capital, data = Grunfeld, type = "honda") 373#' plmtest(inv ~ value + capital, data = Grunfeld, type = "bp") 374#' plmtest(inv ~ value + capital, data = Grunfeld, type = "bp", effect = "twoways") 375#' plmtest(inv ~ value + capital, data = Grunfeld, type = "ghm", effect = "twoways") 376#' plmtest(inv ~ value + capital, data = Grunfeld, type = "kw", effect = "twoways") 377#' 378#' Grunfeld_unbal <- Grunfeld[1:(nrow(Grunfeld)-1), ] # create an unbalanced panel data set 379#' g_unbal <- plm(inv ~ value + capital, data = Grunfeld_unbal, model = "pooling") 380#' plmtest(g_unbal) # unbalanced version of test is indicated in output 381#' 382plmtest <- function(x, ...){ 383 UseMethod("plmtest") 384} 385 386#' @rdname plmtest 387#' @export 388plmtest.plm <- function(x, 389 effect = c("individual", "time", "twoways"), 390 type = c("honda", "bp", "ghm", "kw"), 391 ...) { 392 393 effect <- match.arg(effect) 394 type <- match.arg(type) 395 if (describe(x, "model") != "pooling") x <- update(x, model = "pooling") 396 pdim <- pdim(x) 397 n <- pdim$nT$n 398 T <- pdim$nT$T 399 N_obs <- pdim$nT$N 400 balanced <- pdim$balanced 401 index <- unclass(attr(model.frame(x), "index")) # unclass for speed 402 id <- index[[1L]] 403 time <- index[[2L]] 404 T_i <- pdim$Tint$Ti 405 N_t <- pdim$Tint$nt 406 res <- resid(x) 407 408 ### calc of parts of test statistic ## 409 # calc. is done w/o using matrix calculation, see, e.g., Baltagi/Li (1990), p. 106 410 CP.res <- crossprod(res) 411 A1 <- as.numeric(crossprod(tapply(res, id, sum)) / CP.res - 1) # == A1 <- sum(tapply(res,id,sum)^2) / sum(res^2) - 1 412 A2 <- as.numeric(crossprod(tapply(res, time, sum)) / CP.res - 1) # == A2 <- sum(tapply(res,time,sum)^2) / sum(res^2) - 1 413 414 M11 <- sum(T_i ^ 2) 415 M22 <- sum(N_t ^ 2) 416 417 LM1 <- N_obs * (1 / sqrt(2 * (M11 - N_obs))) * A1 # == sqrt( (((N_obs)^2) / 2) * ( A1^2 / (M11 - N_obs)) ) [except sign due to positive sqrt] 418 LM2 <- N_obs * (1 / sqrt(2 * (M22 - N_obs))) * A2 # == sqrt( (((N_obs)^2) / 2) * ( A2^2 / (M22 - N_obs)) ) [except sign due to positive sqrt] 419 ### END calc of parts of test statistic ## 420 421 422 if (effect != "twoways"){ 423 # oneway 424 if (!type %in% c("honda", "bp", "kw")) 425 stop("type must be one of \"honda\", \"bp\" or \"kw\" for a one way model") # kw oneway coincides with honda 426 427 stat <- if(effect == "individual") LM1 else LM2 428 stat <- switch(type, 429 honda = c(normal = stat), 430 bp = c(chisq = stat ^ 2), 431 kw = c(normal = stat)) 432 433 parameter <- switch(type, 434 honda = NULL, 435 bp = c(df = 1), # df = 1 in the oneway case (Baltagi (2013), p. 70) 436 kw = NULL) 437 438 pval <- switch(type, 439 honda = pnorm(stat, lower.tail = FALSE), # honda oneway ~ N(0,1), alternative is one-sided (Baltagi (2013), p. 71/202) 440 bp = pchisq(stat, df = parameter, lower.tail = FALSE), # df = 1 in the one-way case, alternative is two-sided (Baltagi (2013), p. 70/201) 441 kw = pnorm(stat, lower.tail = FALSE)) # kw oneway ~ N(0,1), alternative is one-sided (Baltagi (2013), p. 71/202) 442 # END oneway 443 } 444 else { # twoways 445 stat <- switch(type, 446 honda = c(normal = (LM1 + LM2) / sqrt(2)), 447 bp = c(chisq = LM1 ^ 2 + LM2 ^ 2), 448 kw = c(normal = (sqrt(M11 - N_obs) / sqrt(M11 + M22 - 2 * N_obs)) * LM1 + 449 (sqrt(M22 - N_obs) / sqrt(M11 + M22 - 2 * N_obs)) * LM2), 450 ghm = c(chibarsq = max(0, LM1) ^ 2 + max(0, LM2) ^ 2)) 451 452 parameter <- switch(type, 453 honda = NULL, 454 bp = c(df = 2), # df = 2 in the twoway case (Baltagi (2013), p. 70/201) 455 kw = NULL, 456 ghm = c(df0 = 0L, df1 = 1L, df2 = 2L, w0 = 1/4, w1 = 1/2, w2 = 1/4)) # chibarsquared (mixed chisq) has several dfs and weights (Baltagi (2013), p. 72/202) 457 458 pval <- switch(type, 459 honda = pnorm(stat, lower.tail = FALSE), # honda two-ways ~ N(0,1), alternative is one-sided (Baltagi (2013), p. 71/202) 460 bp = pchisq(stat, df = parameter, lower.tail = FALSE), # is df = 2 in the twoway case, alternative is two-sided (Baltagi (2013), p. 70/201) 461 kw = pnorm(stat, lower.tail = FALSE), # kw twoways ~ N(0,1), alternative is one-sided (Baltagi (2013), p. 71/202) 462 ghm = pchibarsq(stat, df = c(0L, 1L, 2L), weights = c(1/4, 1/2, 1/4), lower.tail = FALSE)) # mixed chisq (also called chi-bar-square), see Baltagi (2013), pp. 71-72, 74, 88, 202-203, 209 463 } # END twoways 464 465 method.type <- switch(type, 466 honda = "Honda", 467 bp = "Breusch-Pagan", 468 ghm = "Gourieroux, Holly and Monfort", 469 kw = "King and Wu") 470 471 method.effect <- switch(effect, 472 id = "individual effects", 473 time = "time effects", 474 twoways = "two-ways effects") 475 476 balanced.type <- if(balanced) "balanced" else "unbalanced" 477 478 method <- paste("Lagrange Multiplier Test - ", method.effect, 479 " (", method.type, ") for ", balanced.type, " panels", sep="") 480 481 if (type %in% c("honda", "kw")) { 482 RVAL <- list(statistic = stat, 483 p.value = pval, 484 method = method, 485 data.name = data.name(x)) 486 } 487 else { # bp, ghm 488 RVAL <- list(statistic = stat, 489 p.value = pval, 490 method = method, 491 parameter = parameter, 492 data.name = data.name(x)) 493 } 494 495 RVAL$alternative <- "significant effects" # TODO: maybe distinguish b/w one-sided and two-sided alternatives? 496 # (bp: two-sided alt.; all others: one-sided alt.?) 497 498 class(RVAL) <- "htest" 499 return(RVAL) 500} 501 502#' @rdname plmtest 503#' @export 504plmtest.formula <- function(x, data, ..., 505 effect = c("individual", "time", "twoways"), 506 type = c("honda", "bp", "ghm", "kw")) { 507 508 cl <- match.call(expand.dots = TRUE) 509 cl$model <- "pooling" # plmtest is performed on the pooling model... 510 cl$effect <- NULL # ... and pooling model has no argument effect... 511 cl$type <- NULL # ... and no argument type => see below: pass on args effect and type to plmtest.plm() 512 names(cl)[2L] <- "formula" 513 m <- match(plm.arg, names(cl), 0L) 514 cl <- cl[c(1L, m)] 515 cl[[1L]] <- as.name("plm") 516 plm.model <- eval(cl, parent.frame()) 517 plmtest(plm.model, effect = effect, type = type) # pass on args effect and type to plmtest.plm() 518} 519 520 521#' F Test for Individual and/or Time Effects 522#' 523#' Test of individual and/or time effects based on the comparison of the 524#' `within` and the `pooling` model. 525#' 526#' For the `plm` method, the argument of this function is two `plm` 527#' objects, the first being a within model, the second a pooling 528#' model. The effects tested are either individual, time or twoways, 529#' depending on the effects introduced in the within model. 530#' 531#' @aliases pFtest 532#' @param x an object of class `"plm"` or of class `"formula"`, 533#' @param z an object of class `"plm"`, 534#' @param data a `data.frame`, 535#' @param \dots further arguments. 536#' @return An object of class `"htest"`. 537#' @export 538#' @author Yves Croissant 539#' @seealso [plmtest()] for Lagrange multiplier tests of individuals 540#' and/or time effects. 541#' @keywords htest 542#' @examples 543#' 544#' data("Grunfeld", package="plm") 545#' gp <- plm(inv ~ value + capital, data = Grunfeld, model = "pooling") 546#' gi <- plm(inv ~ value + capital, data = Grunfeld, 547#' effect = "individual", model = "within") 548#' gt <- plm(inv ~ value + capital, data = Grunfeld, 549#' effect = "time", model = "within") 550#' gd <- plm(inv ~ value + capital, data = Grunfeld, 551#' effect = "twoways", model = "within") 552#' pFtest(gi, gp) 553#' pFtest(gt, gp) 554#' pFtest(gd, gp) 555#' pFtest(inv ~ value + capital, data = Grunfeld, effect = "twoways") 556#' 557pFtest <- function(x, ...){ 558 UseMethod("pFtest") 559} 560 561#' @rdname pFtest 562#' @export 563pFtest.formula <- function(x, data, ...){ 564 cl <- match.call(expand.dots = TRUE) 565 cl$model <- "within" 566 names(cl)[2L] <- "formula" 567 m <- match(plm.arg,names(cl), 0L) 568 cl <- cl[c(1L, m)] 569 cl[[1L]] <- as.name("plm") 570 plm.within <- eval(cl,parent.frame()) 571 plm.pooling <- update(plm.within, model = "pooling") 572 pFtest(plm.within, plm.pooling, ...) 573} 574 575#' @rdname pFtest 576#' @export 577pFtest.plm <- function(x, z, ...){ 578 within <- x 579 pooling <- z 580 ## leave this interface check commented because pkg AER (reverse dependency) has examples that 581 ## use pFtest(within_twoway, within_time) 582 # if (! (describe(x, "model") == "within" && describe(z, "model") == "pooling")) 583 # stop("the two arguments should be a 'within' and a 'pooling' model (in this order)") 584 585 effect <- describe(x, "effect") 586 df1 <- df.residual(pooling)-df.residual(within) 587 df2 <- df.residual(within) 588 ssrp <- as.numeric(crossprod(residuals(pooling))) 589 ssrw <- as.numeric(crossprod(residuals(within))) 590 stat <- (ssrp-ssrw)/ssrw/df1*df2 591 names(stat) <- "F" 592 parameter <- c(df1, df2) 593 names(parameter) <- c("df1", "df2") 594 pval <- pf(stat, df1, df2, lower.tail = FALSE) 595 alternative <- "significant effects" 596 res <- list(statistic = stat, 597 p.value = pval, 598 method = paste("F test for ", effect, " effects", sep=""), 599 parameter = parameter, 600 data.name = data.name(x), 601 alternative = alternative) 602 class(res) <- "htest" 603 res 604} 605 606############## pwaldtest() ############################################ 607# pwaldtest is used in summary.plm, summary.pht, summary.pgmm to compute the 608# Chi-square or F statistic, but can be used as a stand-alone test of 609# joint significance of all slopes 610# 611# Short intro (but see associated help file) 612# arg 'vcov' non-NULL => the robust tests are carried out 613# arg df2adj == TRUE does finite-sample/cluster adjustment for F tests's df2 614# args .df1, .df2 are only there if user wants to do overwriting of dfs (user has final say) 615# 616# Chi-sq test for IV models as in Wooldridge (1990), A note on the Lagrange multiplier and F-statistics for two stage least 617# squares regressions, Economics Letters 34: 151-155. 618 619#' Wald-style Chi-square Test and F Test 620#' 621#' Wald-style Chi-square test and F test of slope coefficients being 622#' zero jointly, including robust versions of the tests. 623#' 624#' 625#' `pwaldtest` can be used stand--alone with a plm object, a pvcm object, 626#' and a pgmm object (for pvcm objects only the 'random' type is valid and no 627#' further arguments are processed; for pgmm objects only arguments `param` 628#' and `vcov` are valid). It is also used in 629#' [summary.plm()] to produce the F statistic and the Chi-square 630#' statistic for the joint test of coefficients and in [summary.pgmm()]. 631#' 632#' `pwaldtest` performs the test if the slope coefficients of a panel 633#' regression are jointly zero. It does not perform general purpose 634#' Wald-style tests (for those, see [lmtest::waldtest()] (from package 635#' \CRANpkg{lmtest}) or [car::linearHypothesis()] (from package 636#' \CRANpkg{car})). 637#' 638#' If a user specified variance-covariance matrix/function is given in 639#' argument `vcov`, the robust version of the tests are carried out. 640#' In that case, if the F test is requested (`test = "F"`) and no 641#' overwriting of the second degrees of freedom parameter is given (by 642#' supplying argument (`.df2`)), the adjustment of the second degrees 643#' of freedom parameter is performed by default. The second degrees of 644#' freedom parameter is adjusted to be the number of unique elements 645#' of the cluster variable - 1, e. g., the number of individuals minus 1. 646#' For the degrees of freedom adjustment of the F test in general, 647#' see e. g. \insertCite{CAME:MILL:15;textual}{plm}, section VII; 648#' \insertCite{ANDR:GOLS:SCMI:13}{plm}, pp. 126, footnote 4. 649#' 650#' The degrees of freedom adjustment requires the vcov object supplied 651#' or created by a supplied function to carry an attribute called 652#' "cluster" with a known clustering described as a character (for now 653#' this could be either `"group"` or `"time"`). The vcovXX functions 654#' of the package \pkg{plm} provide such an attribute for their 655#' returned variance--covariance matrices. No adjustment is done for 656#' unknown descriptions given in the attribute "cluster" or when the 657#' attribute "cluster" is not present. Robust vcov objects/functions 658#' from package \CRANpkg{clubSandwich} work as inputs to `pwaldtest`'s 659#' F test because a they are translated internally to match the needs 660#' described above. 661#' 662#' @aliases pwaldtest 663#' @param x an estimated model of which the coefficients should be 664#' tested (usually of class `"plm"`/`"pvcm"`/`"pgmm"`)`, 665#' @param test a character, indicating the test to be performed, may 666#' be either `"Chisq"` or `"F"` for the Wald-style 667#' Chi-square test or F test, respectively, 668#' @param vcov `NULL` by default; a `matrix` giving a 669#' variance--covariance matrix or a function which computes such; 670#' if supplied (non `NULL`), the test is carried out using 671#' the variance--covariance matrix indicated resulting in a robust 672#' test, 673#' @param df2adj logical, only relevant for `test = "F"`, 674#' indicating whether the adjustment for clustered standard errors 675#' for the second degrees of freedom parameter should be performed 676#' (see **Details**, also for further requirements regarding 677#' the variance--covariance matrix in `vcov` for the 678#' adjustment to be performed), 679#' @param .df1 a numeric, used if one wants to overwrite the first 680#' degrees of freedom parameter in the performed test (usually not 681#' used), 682#' @param .df2 a numeric, used if one wants to overwrite the second 683#' degrees of freedom parameter for the F test (usually not used), 684#' @param param (for pgmm method only): select the parameters to be tested: 685#' `"coef"`, `"time"`, or `"all"``. 686#' @param \dots further arguments (currently none). 687#' @return An object of class `"htest"`, except for pvcm's within model for which 688#' a data.frame with results of the Wald chi-square tests and F tests per 689#' regression is returned. 690#' @export 691#' @author Yves Croissant (initial implementation) and Kevin Tappe 692#' (extensions: vcov argument and F test's df2 adjustment) 693#' @seealso 694#' 695#' [vcovHC()] for an example of the vcovXX functions, a robust 696#' estimation for the variance--covariance matrix; [summary.plm()] 697#' @references 698#' 699#' \insertRef{WOOL:10}{plm} 700#' 701#' \insertRef{ANDR:GOLS:SCMI:13}{plm} 702#' 703#' \insertRef{CAME:MILL:15}{plm} 704#' 705#' @keywords htest 706#' @examples 707#' 708#' data("Grunfeld", package = "plm") 709#' mod_fe <- plm(inv ~ value + capital, data = Grunfeld, model = "within") 710#' mod_re <- plm(inv ~ value + capital, data = Grunfeld, model = "random") 711#' pwaldtest(mod_fe, test = "F") 712#' pwaldtest(mod_re, test = "Chisq") 713#' 714#' # with robust vcov (matrix, function) 715#' pwaldtest(mod_fe, vcov = vcovHC(mod_fe)) 716#' pwaldtest(mod_fe, vcov = function(x) vcovHC(x, type = "HC3")) 717#' 718#' pwaldtest(mod_fe, vcov = vcovHC(mod_fe), df2adj = FALSE) # w/o df2 adjustment 719#' 720#' # example without attribute "cluster" in the vcov 721#' vcov_mat <- vcovHC(mod_fe) 722#' attr(vcov_mat, "cluster") <- NULL # remove attribute 723#' pwaldtest(mod_fe, vcov = vcov_mat) # no df2 adjustment performed 724#' 725#' 726pwaldtest <- function(x, ...) { 727 UseMethod("pwaldtest") 728} 729 730#' @rdname pwaldtest 731#' @export 732pwaldtest.plm <- function(x, test = c("Chisq", "F"), vcov = NULL, 733 df2adj = (test == "F" && !is.null(vcov) && missing(.df2)), .df1, .df2, ...) { 734 model <- describe(x, "model") 735 test <- match.arg(test) 736 df1 <- if(model == "within") length(coef(x)) else { length(coef(x)) - has.intercept(x) } 737 df2 <- df.residual(x) 738# tss <- tss(x) # not good for models without intercept 739# ssr <- deviance(x) # -- " -- 740 vcov_arg <- vcov 741 int <- "(Intercept)" 742 coefs_wo_int <- coef(x)[!(names(coef(x)) %in% int)] 743 if(!length(coefs_wo_int)) stop(paste("No non-intercept regressors in input model 'x',", 744 "cannot perform Wald joint significance test")) 745 # sanity check 746 if (df2adj == TRUE && (is.null(vcov_arg) || test != "F")) { 747 stop("df2adj == TRUE sensible only for robust F test, i.e., test == \"F\" and !is.null(vcov) and missing(.df2)") 748 } 749 750 # if robust test: prepare robust vcov 751 if (!is.null(vcov_arg)) { 752 if (is.matrix(vcov_arg)) rvcov <- rvcov_orig <- vcov_arg 753 if (is.function(vcov_arg)) rvcov <- rvcov_orig <- vcov_arg(x) 754 755 rvcov_name <- paste0(", vcov: ", paste0(deparse(substitute(vcov)))) # save "name" for later 756 757 if (int %in% names(coef(x))) { # drop intercept, if present 758 rvcov <- rvcov_orig[!rownames(rvcov_orig) %in% int, !colnames(rvcov_orig) %in% int] 759 attr(rvcov, which = "cluster") <- attr(rvcov_orig, which = "cluster") # restore dropped 'cluster' attribute 760 } 761 # if robust F test: by default, do finite-sample adjustment for df2 762 if (df2adj == TRUE && test == "F") { 763 # determine the variable that the clustering is done on by 764 # attribute "cluster" in the vcov (matrix object) 765 # if only one member in cluster: fall back to original df2 766 if (!is.null(attr(rvcov, which = "cluster"))) { 767 768 # if supplied vcov is from package "clubSandwich": translate attr "cluster" to fit our code 769 # (use rvcov_orig here for the test as the above dropping of the intercept drops the special classes of rvcov) 770 if (inherits(rvcov_orig, "vcovCR")) rvcov <- trans_clubSandwich_vcov(CSvcov = rvcov, index = attr(model.frame(x), "index")) 771 772 cluster <- attr(rvcov, which = "cluster") 773 pdim <- pdim(x) 774 df2 <- switch(cluster, 775 group = { if(pdim$nT$n == 1L) df2 else (pdim$nT$n - 1L) }, 776 time = { if(pdim$nT$T == 1L) df2 else (pdim$nT$T - 1L) }, 777 # TODO: what about double clustering? vcovDC? vcovDC identifies itself as attr(obj, "cluster")="group-time") 778 # default: 779 { # warning("unknown/not implemented clustering, no df2 adjustment for finite-samples") 780 df2} 781 ) 782 } else { 783 # no information on clustering found, do not adjust df2 784 # (other options would be: assume cluster = "group", or fall-back to non robust statistics (set vcov_arg <- NULL)) 785 warning("no attribute 'cluster' in robust vcov found, no finite-sample adjustment for df2") # assuming cluster = \"group\"") 786 # df2 <- as.integer(pdim(x)$nT$n - 1) # assume cluster = "group" 787 } 788 } 789 } 790 791 # final say: overwrite Dfs if especially supplied 792 if (!missing(.df1)) df1 <- .df1 793 if (!missing(.df2)) df2 <- .df2 794 795 if (test == "Chisq"){ 796 # perform non-robust chisq test 797 if (is.null(vcov_arg)) { 798 stat <- as.numeric(crossprod(solve(vcov(x)[names(coefs_wo_int), names(coefs_wo_int)], coefs_wo_int), coefs_wo_int)) 799# stat < - (tss-ssr)/(ssr/df2) # does not produce correct results for unbalanced RE models and (un)balanced IV models 800 names(stat) <- "Chisq" 801 pval <- pchisq(stat, df = df1, lower.tail = FALSE) 802 parameter <- c(df = df1) 803 method <- "Wald test for joint significance" 804 } else { 805 # perform robust chisq test 806 stat <- as.numeric(crossprod(solve(rvcov, coefs_wo_int), coefs_wo_int)) 807 names(stat) <- "Chisq" 808 pval <- pchisq(stat, df = df1, lower.tail = FALSE) 809 parameter <- c(df = df1) 810 method <- paste0("Wald test for joint significance (robust)", rvcov_name) 811 } 812 } 813 if (test == "F"){ 814 if(length(formula(x))[2L] > 1L) stop("test = \"F\" not sensible for IV models") 815 if (is.null(vcov_arg)) { 816 # perform "normal" F test 817 stat <- as.numeric(crossprod(solve(vcov(x)[names(coefs_wo_int), names(coefs_wo_int)], coefs_wo_int), coefs_wo_int)) / df1 818# stat <- (tss-ssr)/ssr*df2/df1 # does not produce correct results for unbalanced RE models 819 names(stat) <- "F" 820 pval <- pf(stat, df1 = df1, df2 = df2, lower.tail = FALSE) 821 parameter <- c(df1 = df1, df2 = df2) 822 method <- "F test for joint significance" 823 } else { 824 # perform robust F test 825 stat <- as.numeric(crossprod(solve(rvcov, coefs_wo_int), coefs_wo_int) / df1) 826 names(stat) <- "F" 827 pval <- pf(stat, df1 = df1, df2 = df2, lower.tail = FALSE) 828 parameter <- c(df1 = df1, df2 = df2) 829 method <- paste0("F test for joint significance (robust)", rvcov_name) 830 } 831 } 832 res <- list(data.name = data.name(x), 833 statistic = stat, 834 parameter = parameter, 835 p.value = pval, 836 method = method, 837 alternative = "at least one coefficient is not null" 838 ) 839 class(res) <- "htest" 840 return(res) 841} 842 843#' @rdname pwaldtest 844#' @export 845pwaldtest.pvcm <- function(x, ...) { 846 model <- describe(x, "model") 847 effect <- describe(x, "effect") 848 849 coefs.no.int <- !names(x$coefficients) %in% "(Intercept)" # logical with non-intercept regressors set to TRUE 850 if(!length(names(x$coefficients)[coefs.no.int])) { 851 # error informatively if only-intercept model (no other regressors) 852 stop(paste("No non-intercept regressors in model(s) of input 'x',", 853 "cannot perform Wald joint significance test(s)")) 854 } 855 856 if(model == "within") { 857 # for the within case, simply return a data.frame with all test results 858 # of single estimations (per individual or per time period) 859 860 ii <- switch(effect, "individual" = 1L, "time" = 2L) 861 residl <- split(x$residuals, unclass(index(x))[[ii]]) 862 863 # vcovs and coefficients w/o intercept 864 vcovl <- lapply(x$vcov, function(x) x[coefs.no.int, coefs.no.int]) 865 coefl <- as.list(data.frame(t(x$coefficients[ , coefs.no.int, drop = FALSE]))) 866 df1 <- ncol(x$coefficients[ , coefs.no.int, drop = FALSE]) # ncol is same df1 for all models (as all models estimate the same coefs) 867 df2 <- lengths(residl) - ncol(x$coefficients) # (any intercept is subtracted) 868 869 statChisqs <- mapply(FUN = function(v, c) as.numeric(crossprod(solve(v, c), c)), 870 vcovl, coefl) 871 statFs <- statChisqs / df1 872 873 pstatChisqs <- pchisq(statChisqs, df = df1, lower.tail = FALSE) 874 pstatFs <- pf(statFs, df1 = df1, df2 = df2, lower.tail = FALSE) 875 876 stats.pvcm.within <- as.data.frame(cbind("Chisq" = statChisqs, 877 "p(chisq)" = pstatChisqs, 878 "F" = statFs, 879 "p(F)" = pstatFs, 880 "df1" = rep(df1, length(residl)), 881 "df2" = df2)) 882 # early return 883 return(stats.pvcm.within) 884 } 885 886 ## case: model == "random" 887 coefs_wo_int <- x$coefficients[coefs.no.int] 888 stat <- as.numeric(crossprod(solve(vcov(x)[coefs.no.int, coefs.no.int], coefs_wo_int), coefs_wo_int)) 889 names(stat) <- "Chisq" 890 df1 <- length(coefs_wo_int) 891 pval <- pchisq(stat, df = df1, lower.tail = FALSE) 892 parameter <- c(df = df1) 893 method <- "Wald test for joint significance" 894 895 res <- list(data.name = data.name(x), 896 statistic = stat, 897 parameter = parameter, 898 p.value = pval, 899 method = method, 900 alternative = "at least one coefficient is not null" 901 ) 902 class(res) <- "htest" 903 return(res) 904} 905 906 907#' @rdname pwaldtest 908#' @export 909pwaldtest.pgmm <- function(x, param = c("coef", "time", "all"), vcov = NULL, ...) { 910 param <- match.arg(param) 911 vcov_supplied <- !is.null(vcov) 912 myvcov <- vcov 913 if (is.null(vcov)) vv <- vcov(x) 914 else if (is.function(vcov)) vv <- myvcov(x) 915 else vv <- myvcov 916 917 model <- describe(x, "model") 918 effect <- describe(x, "effect") 919 if (param == "time" && effect == "individual") stop("no time dummies in this model") 920 transformation <- describe(x, "transformation") 921 if (model == "onestep") coefficients <- x$coefficients 922 else coefficients <- x$coefficients[[2L]] 923 Ktot <- length(coefficients) 924 Kt <- length(x$args$namest) 925 926 switch(param, 927 "time" = { 928 start <- Ktot - Kt + if(transformation == "ld") 2 else 1 929 end <- Ktot 930 }, 931 "coef" = { 932 start <- 1 933 end <- if (effect == "twoways") Ktot - Kt else Ktot 934 }, 935 "all" = { 936 start <- 1 937 end <- Ktot 938 }) 939 coef <- coefficients[start:end] 940 vv <- vv[start:end, start:end] 941 stat <- as.numeric(crossprod(coef, crossprod(solve(vv), coef))) 942 names(stat) <- "chisq" 943 parameter <- length(coef) 944 names(parameter) <- "df" 945 pval <- pchisq(stat, df = parameter, lower.tail = FALSE) 946 method <- "Wald test for joint significance" 947 if (vcov_supplied) { 948 rvcov_name <- paste0(", vcov: ", paste0(deparse(substitute(vcov)))) 949 method <- paste0(method, " (robust)", rvcov_name) 950 } 951 wald.pgmm <- list(statistic = stat, 952 p.value = pval, 953 parameter = parameter, 954 method = method, 955 alternative = "at least one coefficient is not null", 956 data.name = data.name(x)) 957 class(wald.pgmm) <- "htest" 958 return(wald.pgmm) 959} 960 961pwaldtest.default <- function(x, ...) { 962 pwaldtest.plm(x, ...) 963} 964 965 966# trans_clubSandwich_vcov: helper function for pwaldtest() 967# translate vcov object from package clubSandwich so it is suitable for summary.plm, plm's pwaldtest. 968# Attribute "cluster" in clubSandwich's vcov objects contains the cluster variable itself. 969# plm's vcov object also has attribute "cluster" but it contains a character as 970# information about the cluster dimension (either "group" or "time") 971# 972# inputs: 973# * CSvcov: a vcov as returned by clubSandwich's vcovCR function [class c("vcovCR", "clubSandwich")] 974# * index: the index belonging to a plm object/model 975# return value: 976# * modified CSvcov (substituted attribute "cluster" with suitable character or NULL) 977trans_clubSandwich_vcov <- function(CSvcov, index) { 978 clustervar <- attr(CSvcov, "cluster") 979 if (!is.null(clustervar)) { 980 if (isTRUE(all.equal(index[[1L]], clustervar))) { 981 attr(CSvcov, "cluster") <- "group" 982 return(CSvcov) 983 } 984 if (isTRUE(all.equal(index[[2L]], clustervar))) { 985 attr(CSvcov, "cluster") <- "time" 986 return(CSvcov) 987 } else { 988 attr(CSvcov, "cluster") <- NULL 989 return(CSvcov) 990 } 991 } 992 warning("no attribute \"cluster\" found in supplied vcov object") 993 return(CSvcov) 994} 995 996 997 998#' Test of Poolability 999#' 1000#' A Chow test for the poolability of the data. 1001#' 1002#' `pooltest` is a *F* test of stability (or Chow test) for the 1003#' coefficients of a panel model. For argument `x`, the estimated 1004#' `plm` object should be a `"pooling"` model or a `"within"` model 1005#' (the default); intercepts are assumed to be identical in the first 1006#' case and different in the second case. 1007#' 1008#' @aliases pooltest 1009#' @param x an object of class `"plm"` for the plm method; an object of 1010#' class `"formula"` for the formula interface, 1011#' @param z an object of class `"pvcm"` obtained with 1012#' `model="within"`, 1013#' @param data a `data.frame`, 1014#' @param \dots further arguments passed to plm. 1015#' @return An object of class `"htest"`. 1016#' @export 1017#' @author Yves Croissant 1018#' @keywords htest 1019#' @examples 1020#' 1021#' data("Gasoline", package = "plm") 1022#' form <- lgaspcar ~ lincomep + lrpmg + lcarpcap 1023#' gasw <- plm(form, data = Gasoline, model = "within") 1024#' gasp <- plm(form, data = Gasoline, model = "pooling") 1025#' gasnp <- pvcm(form, data = Gasoline, model = "within") 1026#' pooltest(gasw, gasnp) 1027#' pooltest(gasp, gasnp) 1028#' 1029#' pooltest(form, data = Gasoline, effect = "individual", model = "within") 1030#' pooltest(form, data = Gasoline, effect = "individual", model = "pooling") 1031#' 1032pooltest <- function(x,...){ 1033 UseMethod("pooltest") 1034} 1035 1036 1037#' @rdname pooltest 1038#' @export 1039pooltest.plm <- function(x, z, ...){ 1040 rss <- deviance(x) 1041 uss <- as.numeric(crossprod(residuals(z))) 1042 dlr <- df.residual(x) 1043 dlu <- df.residual(z) 1044 df1 <- dlr - dlu 1045 df2 <- dlu 1046 stat <- (rss-uss)/uss*df2/df1 1047 pval <- pf(stat, df1 = df1, df2 = df2, lower.tail = FALSE) 1048 parameter <- c(df1 = df1, df2 = df2) 1049 names(stat) <- "F" 1050 res <- list(statistic = stat, 1051 parameter = parameter, 1052 p.value = pval, 1053 data.name = data.name(x), 1054 alternative = "unstability", 1055 method = "F statistic") 1056 class(res) <- "htest" 1057 res 1058} 1059 1060#' @rdname pooltest 1061#' @export 1062pooltest.formula <- function(x, data, ...){ 1063 cl <- match.call(expand.dots = TRUE) 1064 cl[[1L]] <- as.name("plm") 1065 names(cl)[[2L]] <- "formula" 1066 if (is.null(cl$effect)) cl$effect <- "individual" 1067 plm.model <- eval(cl, parent.frame()) 1068 1069 cl[[1L]] <- as.name("pvcm") 1070 names(cl)[[2L]] <- "formula" 1071 if (is.null(cl$effect)) cl$effect <- "individual" 1072 cl$model <- "within" 1073 pvcm.model <- eval(cl, parent.frame()) 1074 1075 pooltest(plm.model, pvcm.model) 1076} 1077 1078