1 2############## Pesaran's CD test and Breusch/Pagan LM Test (also scaled) ############### 3 4 ## Pesaran's CD test for cross-sectional dependence in panel data models 5 ## (and Breusch and Pagan's LM and scaled LM) 6 ## ref. Pesaran, General diagnostic tests..., CESifo WP 1229, 2004 7 8 ## In case K+1>T the group-specific model is not estimable; 9 ## as in Greene 11.7.2, formula (11.23) we use the group-specific residuals 10 ## of a consistent estimator. This may be pooled OLS, RE, FE. Here the 11 ## default is set to FE. 12 13 ## Note that the test can be performed on the results of plm objects with 14 ## any kind of effects: having "time" effects means checking for 15 ## xs-dependence *after* introducing time dummies. 16 17 ## In principle, the test can be performed on the results of *any* 18 ## panelmodel object. Some issues remain regarding standardization of 19 ## model output: some missing pieces are, e.g., the 'model$indexes' 20 ## in ggls. ''fd'' models are also not compatible because of indexes 21 ## keeping the original timespan, while data lose the first period. 22 23## production version, generic and based on plm 24 25## version 11: added test = "bcsclm" 26## 27## version 10: 28## substantial optimization for speed, now fast (few seconds) on N=3000 29## all methods pass on a pseries to pcdres() 30 31## make toy example 32#dati <- data.frame(ind=rep(1:7, 4), time=rep(1:4, each=7), x=rnorm(28), 33# group=rep(c(1,1,2,2,2,3,3), 4)) 34#pdati <- pdata.frame(dati) 35 36#' Tests of cross-section dependence for panel models 37#' 38#' Pesaran's CD or Breusch--Pagan's LM (local or global) tests for cross 39#' sectional dependence in panel models 40#' 41#' These tests are originally meant to use the residuals of separate 42#' estimation of one time--series regression for each cross-sectional 43#' unit in order to check for cross--sectional dependence (`model = NULL`). 44#' If a different model specification (`model = "within"`, `"random"`, 45#' \ldots{}) is assumed consistent, one can resort to its residuals for 46#' testing (which is common, e.g., when the time dimension's length is 47#' insufficient for estimating the heterogeneous model). 48#' 49#' If the time 50#' dimension is insufficient and `model = NULL`, the function defaults 51#' to estimation of a `within` model and issues a warning. The main 52#' argument of this function may be either a model of class 53#' `panelmodel` or a `formula` and `data frame`; in the second case, 54#' unless `model` is set to `NULL`, all usual parameters relative to 55#' the estimation of a `plm` model may be passed on. The test is 56#' compatible with any consistent `panelmodel` for the data at hand, 57#' with any specification of `effect` (except for `test = "bcsclm"` which 58#' requires a within model with either individual or two-ways effect). 59#' E.g., specifying `effect = "time"` or `effect = "twoways"` allows 60#' to test for residual cross-sectional dependence after the introduction 61#' of time fixed effects to account for common shocks. 62#' 63#' A **local** version of either test can be computed by supplying a 64#' proximity matrix (elements coercible to `logical`) with argument 65#' `w` which provides information on whether any pair of individuals 66#' are neighbours or not. If `w` is supplied, only neighbouring pairs 67#' will be used in computing the test; else, `w` will default to 68#' `NULL` and all observations will be used. The matrix need not be 69#' binary, so commonly used "row--standardized" matrices can be 70#' employed as well. `nb` objects from \CRANpkg{spdep} must instead be 71#' transformed into matrices by \CRANpkg{spdep}'s function `nb2mat` 72#' before using. 73#' 74#' The methods implemented are suitable also for unbalanced panels. 75#' 76#' Pesaran's CD test (`test="cd"`), Breusch and Pagan's LM test 77#' (`test="lm"`), and its scaled version (`test="sclm"`) are all 78#' described in \insertCite{PESA:04;textual}{plm} (and complemented by 79#' Pesaran (2005)). The bias-corrected scaled test (`test="bcsclm"`) 80#' is due to \insertCite{BALT:FENG:KAO:12}{plm} and only valid for 81#' within models including the individual effect (it's unbalanced 82#' version uses max(Tij) for T) in the bias-correction term). 83#' \insertCite{BREU:PAGA:80;textual}{plm} is the original source for 84#' the LM test. 85#' 86#' The test on a `pseries` is the same as a test on a pooled 87#' regression model of that variable on a constant, i.e., 88#' `pcdtest(some_pseries)` is equivalent to `pcdtest(plm(some_var ~ 1, 89#' data = some_pdata.frame, model = "pooling")` and also equivalent to 90#' `pcdtest(some_var ~ 1, data = some_data)`, where `some_var` is 91#' the variable name in the data which corresponds to `some_pseries`. 92#' 93#' @aliases pcdtest 94#' @param x an object of class `formula`, `panelmodel`, or `pseries` 95#' (depending on the respective interface) describing the model to 96#' be tested, 97#' @param data a `data.frame`, 98#' @param index an optional numerical index, if `NULL`, the first two 99#' columns of the data.frame provided in argument `data` are 100#' assumed to be the index variables; for further details see 101#' [pdata.frame()], 102#' @param model an optional character string indicating which type of 103#' model to estimate; if left to `NULL`, the original 104#' heterogeneous specification of Pesaran is used, 105#' @param test the type of test statistic to be returned. One of 106#' \itemize{ \item `"cd"` for Pesaran's CD statistic, \item `"lm"` 107#' for Breusch and Pagan's original LM statistic, \item `"sclm"` 108#' for the scaled version of Breusch and Pagan's LM statistic, 109#' \item `"bcsclm"` for the bias-corrected scaled version of 110#' Breusch and Pagan's LM statistic, \item `"rho"` for the average 111#' correlation coefficient, \item `"absrho"` for the average 112#' absolute correlation coefficient,} 113#' @param w either `NULL` (default) for the global tests or -- for the 114#' local versions of the statistics -- a `n x n` `matrix` 115#' describing proximity between individuals, with \eqn{w_ij = a} 116#' where \eqn{a} is any number such that `as.logical(a)==TRUE`, if 117#' \eqn{i,j} are neighbours, \eqn{0} or any number \eqn{b} such 118#' that `as.logical(b)==FALSE` elsewhere. Only the lower 119#' triangular part (without diagonal) of `w` after coercing by 120#' `as.logical()` is evaluated for neighbouring information (but 121#' `w` can be symmetric). See also **Details** and 122#' **Examples**, 123#' @param \dots further arguments to be passed on for model estimation to `plm`, 124#' such as `effect` or `random.method`. 125#' @return An object of class `"htest"`. 126#' @export 127#' @references 128#' 129#' \insertRef{BALT:FENG:KAO:12}{plm} 130#' 131#' \insertRef{BREU:PAGA:80}{plm} 132#' 133#' \insertRef{PESA:04}{plm} 134#' 135#' \insertRef{PESA:15}{plm} 136#' 137#' @keywords htest 138#' @examples 139#' 140#' data("Grunfeld", package = "plm") 141#' ## test on heterogeneous model (separate time series regressions) 142#' pcdtest(inv ~ value + capital, data = Grunfeld, 143#' index = c("firm", "year")) 144#' 145#' ## test on two-way fixed effects homogeneous model 146#' pcdtest(inv ~ value + capital, data = Grunfeld, model = "within", 147#' effect = "twoways", index = c("firm", "year")) 148#' 149#' ## test on panelmodel object 150#' g <- plm(inv ~ value + capital, data = Grunfeld, index = c("firm", "year")) 151#' pcdtest(g) 152#' 153#' ## scaled LM test 154#' pcdtest(g, test = "sclm") 155#' 156#' ## test on pseries 157#' pGrunfeld <- pdata.frame(Grunfeld) 158#' pcdtest(pGrunfeld$value) 159#' 160#' ## local test 161#' ## define neighbours for individual 2: 1, 3, 4, 5 in lower triangular matrix 162#' w <- matrix(0, ncol= 10, nrow=10) 163#' w[2,1] <- w[3,2] <- w[4,2] <- w[5,2] <- 1 164#' pcdtest(g, w = w) 165#' 166pcdtest <- function(x, ...) 167{ 168 UseMethod("pcdtest") 169} 170 171## this formula method here only for adding "rho" and "absrho" 172## arguments 173 174#' @rdname pcdtest 175#' @export 176pcdtest.formula <- function(x, data, index = NULL, model = NULL, 177 test = c("cd", "sclm", "bcsclm", "lm", "rho", "absrho"), 178 w = NULL, ...) { 179 #data <- pdata.frame(data, index = index) 180 test <- match.arg(test) 181 if(test == "bcsclm" && (is.null(model) || model != "within")) 182 stop("for test = 'bcsclm', set argument model = 'within'") 183 184 # evaluate formula in parent frame 185 cl <- match.call(expand.dots = TRUE) 186 cl$model <- if(test != "bcsclm") "pooling" else "within" 187 if(test == "bcsclm") { 188 # check args model and effect for test = "bcsclm" 189 if(is.null(cl$effect)) cl$effect <- "individual" # make default within model is individual within 190 eff <- isTRUE(cl$effect == "individual" || cl$effect == "twoways") 191 if(model != "within" || !eff) stop("for test = 'bcsclm', requirement is model = \"within\" and effect = \"individual\" or \"twoways\"") 192 } 193 names(cl)[2L] <- "formula" 194 m <- match(plm.arg, names(cl), 0L) 195 cl <- cl[c(1L, m)] 196 cl[[1L]] <- as.name("plm") 197 mymod <- eval(cl, parent.frame()) # mymod is either "pooling" or "within" (the latter iff for test = "bcsclm") 198 199 hetero.spec <- if(is.null(model)) TRUE else FALSE 200 201 if(hetero.spec && min(pdim(mymod)$Tint$Ti) < length(mymod$coefficients)+1) { 202 warning("Insufficient number of observations in time to estimate heterogeneous model: using within residuals", 203 call. = FALSE) 204 hetero.spec <- FALSE 205 model <- "within" 206 } 207 208 ind0 <- attr(model.frame(mymod), "index") 209 tind <- as.numeric(ind0[[2L]]) 210 ind <- as.numeric(ind0[[1L]]) 211 212 if(hetero.spec) { 213 ## estimate individual normal regressions one by one 214 ## (original heterogeneous specification of Pesaran) 215 X <- model.matrix(mymod) 216 y <- model.response(model.frame(mymod)) 217 unind <- unique(ind) 218 n <- length(unind) 219 ti.res <- vector("list", n) 220 ind.res <- vector("list", n) 221 tind.res <- vector("list", n) 222 for (i in 1:n) { 223 tX <- X[ind == unind[i], , drop = FALSE] 224 ty <- y[ind == unind[i]] 225 res.i <- lm.fit(tX, ty)$residuals 226 ti.res[[i]] <- res.i 227 names(ti.res[[i]]) <- tind[ind == unind[i]] 228 ind.res[[i]] <- rep(i, length(res.i)) 229 tind.res[[i]] <- tind[ind == unind[i]] 230 } 231 ## make pseries of (all) residuals 232 resdata <- data.frame(ee = unlist(ti.res, use.names = FALSE), 233 ind = unlist(ind.res, use.names = FALSE), 234 tind = unlist(tind.res, use.names = FALSE)) 235 pee <- pdata.frame(resdata, index = c("ind", "tind")) 236 tres <- pee$ee 237 } else { 238 # else case is one of: 239 # a) insufficient number of observations for heterogen. spec. or 240 # b) model specified when function was called (incl. case test = "bcsclm") 241 if(test != "bcsclm") { 242 # Estimate the model specified originally in function call or due to 243 # forced model switch to within model by insufficient number of 244 # observations for heterogen. spec. 245 # (for test = "bcsclm" it is ensured that a within model was already 246 # estimated -> no need to estimate again a within model) 247 cl$model <- model 248 mymod <- eval(cl, parent.frame()) 249 } 250 251 tres <- resid(mymod) 252 unind <- unique(ind) 253 n <- length(unind) 254 t <- min(pdim(mymod)$Tint$Ti) 255 nT <- length(ind) 256 k <- length(mymod$coefficients) 257 } 258 259 return(pcdres(tres = tres, n = n, w = w, 260 form = paste(deparse(x)), 261 test = test)) 262} 263 264 265## panelmodel method: just fetch resid (as a pseries) and hand over to pcdres 266 267#' @rdname pcdtest 268#' @export 269pcdtest.panelmodel <- function(x, test = c("cd", "sclm", "bcsclm", "lm", "rho", "absrho"), 270 w = NULL, ...) { 271 272 test <- match.arg(test) 273 model <- describe(x, "model") 274 effect <- describe(x, "effect") 275 eff <- (effect == "individual" || effect == "twoways") 276 if (test == "bcsclm") 277 if (model != "within" || !eff) stop("for test = 'bcsclm', model x must be a within individual or twoways model") 278 279 tres <- resid(x) 280 index <- attr(model.frame(x), "index") 281 #tind <- as.numeric(index[[2L]]) 282 ind <- as.numeric(index[[1L]]) 283 unind <- unique(ind) 284 n <- length(unind) 285 #t <- pdim(x)$Tint$Ti 286 #nT <- length(ind) 287 #k <- length(x$coefficients) 288 return(pcdres(tres = tres, n = n, w = w, 289 form = paste(deparse(x$formula)), 290 test = test)) 291} 292 293#' @rdname pcdtest 294#' @export 295pcdtest.pseries <- function(x, test = c("cd", "sclm", "bcsclm", "lm", "rho", "absrho"), 296 w = NULL, ...) { 297 298 ## calculates local or global CD test on a pseries 'x' just as it 299 ## would on model residuals 300 ## important difference here: a pseries _can_ have NAs 301 302 # input check 303 if (!inherits(x, "pseries")) stop("input 'x' needs to be of class \"pseries\"") 304 form <- paste(deparse(substitute(x))) 305 306 pos.na <- is.na(x) 307 if (any(pos.na)) { 308 x <- subset_pseries(x, !pos.na) # TODO: use [.pseries (pseries subsetting) once implemented 309 warning("NA values encountered in input and removed") 310 if (length(x) == 0L) stop("input is empty after removal of NA values") 311 } 312 313 ## get indices 314 tind <- as.numeric(attr(x, "index")[[2L]]) 315 ind <- as.numeric(attr(x, "index")[[1L]]) 316 317 ## det. number of groups and df 318 unind <- unique(ind) 319 n <- length(unind) 320 321 tres <- x 322 323 ## "pre-allocate" an empty list of length n 324 #tres <- vector("list", n) 325 326 ## use model residuals, group by group 327 ## list of n: 328 ## t_i residuals for each x-sect. 1..n 329 #for(i in 1:n) { 330 # # remove NAs 331 # xnonna <- !is.na(x[ind==unind[i]]) 332 # tres[[i]] <- x[ind==unind[i]][xnonna] 333 # ## name resids after the time index 334 # names(tres[[i]]) <- tind[ind==unind[i]][xnonna] 335 # } 336 337 return(pcdres(tres = tres, n = n, w = w, 338 form = form, 339 test = match.arg(test))) 340} 341 342pcdres <- function(tres, n, w, form, test) { 343 # 'form' is a character describing the formula (not a formula object!) 344 # and goes into htest_object$data.name 345 346 ## Take model residuals as pseries, and calc. test 347 ## (from here on, what's needed for rho_ij is ok) 348 349 ## this function is the modulus calculating the test, 350 ## to be called from pcdtest.formula, 351 ## pcdtest.panelmodel or pcdtest.pseries 352 353 ## now (since v10) tres is the pseries of model residuals 354 355 ## calc matrix of all possible pairwise corr. 356 ## coeffs. (200x speedup from using cor()) 357 wideres <- t(preshape(tres, na.rm = FALSE)) 358 rho <- cor(wideres, use = "pairwise.complete.obs") 359 360 ## find length of intersecting pairs 361 ## fast method, times down 200x 362 data.res <- data.frame(time = attr(tres, "index")[[2L]], 363 indiv = attr(tres, "index")[[1L]]) 364 ## tabulate which obs in time for each ind are !na 365 presence.tab <- table(data.res) 366 ## calculate t.ij 367 t.ij <- crossprod(presence.tab) 368 369 # input check 370 if (!is.null(w)) { 371 dims.w <- dim(w) 372 if(dims.w[1L] != n || dims.w[2L] != n) 373 stop(paste0("matrix 'w' describing proximity of individuals has wrong dimensions: ", 374 "should be ", n, " x ", n, " (no. of individuals) but is ", dims.w[1L], " x ", dims.w[2L])) 375 } 376 377 378 ## begin features for local test #################### 379 ## higher orders omitted for now, use wlag() explicitly 380 381 ## if global test, set all elements in w to 1 382 if(is.null(w)) { 383 w <- matrix(1, ncol = n, nrow = n) 384 dep <- "" 385 } else { dep <- "local" } 386 387 ## make (binary) selector matrix based on the contiguity matrix w 388 ## and extracting elements corresponding to ones in the lower triangle 389 ## excluding the diagonal 390 391 ## transform in logicals (0=FALSE, else=TRUE: no need to worry 392 ## about row-std. matrices) 393 selector.mat <- matrix(as.logical(w), ncol = n) 394 395 ## some sanity checks for 'w' (not perfect sanity, but helps) 396 if (sum(selector.mat[lower.tri(selector.mat, diag = FALSE)]) == 0) { 397 stop(paste0("no neighbouring individuals defined in proximity matrix 'w'; ", 398 "only lower triangular part of 'w' (w/o diagonal) is evaluated")) 399 } else { 400 if (sum(selector.mat[upper.tri(selector.mat, diag = FALSE)]) != 0) { 401 if (!isSymmetric((unname(selector.mat)))) { # unname needed to ignore rownames and colnames 402 stop(paste0("proximity matrix 'w' is ambiguous: upper and lower triangular part ", 403 "define different neighbours (it is sufficient to provide information ", 404 "about neighbours only in the lower triangluar part of 'w'")) 405 } 406 } 407 } 408 409 ## if no intersection or only 1 shared period of e_it and e_jt 410 ## => exclude from calculation and issue a warning. 411 ## In general, length(m.ij) gives the number of shared periods by indiviudals i, j 412 ## Thus, non intersecting pairs are indicated by length(m.ij) == 0 (t.ij[i,j] == 0) 413 no.one.intersect <- (t.ij <= 1) 414 if (any(no.one.intersect, na.rm = TRUE)) { 415 # t.ij is a lower triangular matrix: do not divide by 2 to get the number of non-intersecting pairs! 416 number.of.non.one.intersecting.pairs <- sum(no.one.intersect, na.rm = TRUE) 417 number.of.total.pairs <- (n*(n-1))/2 418 share.on.one.intersect.pairs <- number.of.non.one.intersecting.pairs / number.of.total.pairs * 100 419 warning(paste("Some pairs of individuals (", 420 signif(share.on.one.intersect.pairs, digits = 2), 421 " percent) do not have any or just one time period in common and have been omitted from calculation", sep="")) 422 selector.mat[no.one.intersect] <- FALSE 423 } 424 425 ## set upper tri and diagonal to FALSE 426 selector.mat[upper.tri(selector.mat, diag = TRUE)] <- FALSE 427 428 ## number of elements in selector.mat 429 ## elem.num = 2*(N*(N-1)) in Pesaran (2004), formulae (6), (7), (31), ... 430 elem.num <- sum(selector.mat) 431 432 ## end features for local test ###################### 433 434 ## Breusch-Pagan or Pesaran statistic for cross-sectional dependence, 435 ## robust vs. unbalanced panels: 436 437 switch(test, 438 lm = { 439 CDstat <- sum((t.ij*rho^2)[selector.mat]) 440 pCD <- pchisq(CDstat, df = elem.num, lower.tail = FALSE) 441 names(CDstat) <- "chisq" 442 parm <- elem.num 443 names(parm) <- "df" 444 testname <- "Breusch-Pagan LM test" 445 }, 446 sclm = { 447 CDstat <- sqrt(1/(2*elem.num))*sum((t.ij*rho^2-1)[selector.mat]) 448 pCD <- 2*pnorm(abs(CDstat), lower.tail = FALSE) 449 names(CDstat) <- "z" 450 parm <- NULL 451 testname <- "Scaled LM test" 452 }, 453 bcsclm = { 454 # Baltagi/Feng/Kao (2012), formula (11) 455 # (unbalanced case as sclm + in bias correction as EViews: max(T_ij) instead of T) 456 CDstat <- sqrt(1/(2*elem.num))*sum((t.ij*rho^2-1)[selector.mat]) - (n/(2*(max(t.ij)-1))) 457 pCD <- 2*pnorm(abs(CDstat), lower.tail = FALSE) 458 names(CDstat) <- "z" 459 parm <- NULL 460 testname <- "Bias-corrected Scaled LM test" 461 }, 462 cd = { 463 # (Pesaran (2004), formula (31)) 464 CDstat <- sqrt(1/elem.num)*sum((sqrt(t.ij)*rho)[selector.mat]) 465 pCD <- 2*pnorm(abs(CDstat), lower.tail = FALSE) 466 names(CDstat) <- "z" 467 parm <- NULL 468 testname <- "Pesaran CD test" 469 }, 470 rho = { 471 CDstat <- sum(rho[selector.mat])/elem.num 472 pCD <- NULL 473 names(CDstat) <- "rho" 474 parm <- NULL 475 testname <- "Average correlation coefficient" 476 }, 477 absrho = { 478 CDstat <- sum(abs(rho)[selector.mat])/elem.num 479 pCD <- NULL 480 names(CDstat) <- "|rho|" 481 parm <- NULL 482 testname <- "Average absolute correlation coefficient" 483 }) 484 485 ##(insert usual htest features) 486 RVAL <- list(statistic = CDstat, 487 parameter = parm, 488 method = paste(testname, "for", dep, 489 "cross-sectional dependence in panels"), 490 alternative = "cross-sectional dependence", 491 p.value = pCD, 492 data.name = form) 493 class(RVAL) <- "htest" 494 return(RVAL) 495} 496 497preshape <- function(x, na.rm = TRUE, ...) { 498 ## reshapes pseries, 499 ## e.g., of residuals from a panelmodel, 500 ## in wide form 501 inames <- names(attr(x, "index")) 502 mres <- reshape(cbind(as.vector(x), attr(x, "index")), 503 direction = "wide", 504 timevar = inames[2L], 505 idvar = inames[1L]) 506 ## drop ind in first column 507 mres <- mres[ , -1L, drop = FALSE] 508 ## reorder columns (may be scrambled depending on first 509 ## available obs in unbalanced panels) 510 mres <- mres[ , order(dimnames(mres)[[2L]])] 511 ## if requested, drop columns (time periods) with NAs 512 if(na.rm) { 513 na.cols <- vapply(mres, FUN = anyNA, FUN.VALUE = TRUE, USE.NAMES = FALSE) 514 if(sum(na.cols) > 0L) mres <- mres[ , !na.cols] 515 } 516 return(mres) 517} 518 519 520 521 522#' Cross--sectional correlation matrix 523#' 524#' Computes the cross--sectional correlation matrix 525#' 526#' 527#' @param x an object of class `pseries` 528#' @param grouping grouping variable, 529#' @param groupnames a character vector of group names, 530#' @param value to complete, 531#' @param \dots further arguments. 532#' @return A matrix with average correlation coefficients within a group 533#' (diagonal) and between groups (off-diagonal). 534#' @export 535#' @keywords htest 536#' @examples 537#' 538#' data("Grunfeld", package = "plm") 539#' pGrunfeld <- pdata.frame(Grunfeld) 540#' grp <- c(rep(1, 100), rep(2, 50), rep(3, 50)) # make 3 groups 541#' cortab(pGrunfeld$value, grouping = grp, groupnames = c("A", "B", "C")) 542#' 543cortab <- function(x, grouping, groupnames = NULL, 544 value = "statistic", ...) { 545 ## makes matrix containing within (diagonal) and between (off-diagonal) 546 ## correlation 547 ## needs a pseries and a groupings vector of **same length** 548 549 ## would use a better naming, and also passing a char or factor as 550 ## grouping index 551 552 ## x must be a pseries 553 if(!inherits(x, "pseries")) stop("First argument must be a pseries") 554 if(length(x) != length(grouping)) stop("arguments 'x' and 'grouping' must have same length") 555 556 fullind <- as.numeric(attr(x, "index")[ , 1L]) 557 ids <- unique(fullind) 558 n <- length(ids) 559 regs <- 1:length(unique(grouping)) 560 561 if(!(is.numeric(grouping))) grouping <- as.numeric(as.factor(grouping)) 562 563 idnames <- as.character(ids) 564 if(is.null(groupnames)) { 565 groupnames <- as.character(unique(grouping)) 566 } 567 568 ## make matrices of between-regions correlations 569 ## (includes within correlation on diagonal) 570 ## for each pair of regions (nb: no duplicates, e.g., 3.1 but not 1.3) 571 572 ## make w<1.n>: 573 for(h in 1:length(regs)) { 574 for(k in 1:h) { 575 statew <- matrix(0, ncol = n, nrow = n) 576 ## make statew for cor. between h and k 577 for(i in 1:n) { 578 ## get first region (all values equal, so take first one) 579 ireg <- grouping[fullind == ids[i]][1L] 580 if(ireg == h) { 581 for(j in 1:n) { 582 jreg <- grouping[fullind == ids[j]][1L] 583 if(jreg == k) statew[i, j] <- 1 584 } 585 } 586 } 587 if(h!=k) statew <- statew + t(statew) 588 ## just for debugging reasons: 589 dimnames(statew) <- list(idnames, idnames) 590 ## eliminate self.correlation of states if i=j 591 diag(statew) <- 0 592 ## not needed: pcdtest seems to do this by construction 593 eval(parse(text=paste("w", h, ".", k, " <- statew", sep=""))) 594 } 595 } 596 597 ## notice: without the line 598 ## '' if(i!=j) statew <- statew + t(statew) '' 599 ## all wn.n matrices would have values only on one half (upper 600 ## or lower triangle) 601 602 ## make generic table of regions' within and between correlation 603 ## argument: a pseries 604 #YC regnames is undefined, so is myw 605 tab.g <- function(x, regs, regnames, test="rho", value) { 606 myw <- 0 607 tabg <- matrix(NA, ncol=length(regs), nrow=length(regs)) 608 for(i in 1:length(regs)) { 609 for(j in 1:i) { 610 ## take appropriate w matrix 611 eval(parse(text = paste("myw<-w", i, ".", j, sep = ""))) 612 tabg[i, j] <- pcdtest(x, test = "rho", w = myw)[[value]] 613 } 614 } 615 dimnames(tabg) <- list(groupnames, groupnames) 616 return(tabg) 617 } 618 regnames <- "" 619 mytab <- tab.g(x, regs = regs, regnames = regnames, test = "rho", value = value) 620 return(mytab) 621} 622 623 624 625 626