1#' Panel estimators for limited dependent variables 2#' 3#' Fixed and random effects estimators for truncated or censored 4#' limited dependent variable 5#' 6#' `pldv` computes two kinds of models: a LSQ/LAD estimator for the 7#' first-difference model (`model = "fd"`) and a maximum likelihood estimator 8#' with an assumed normal distribution for the individual effects 9#' (`model = "random"` or `"pooling"`). 10#' 11#' For maximum-likelihood estimations, `pldv` uses internally function 12#' [maxLik::maxLik()] (from package \CRANpkg{maxLik}). 13#' 14#' @aliases pldv 15#' @param formula a symbolic description for the model to be 16#' estimated, 17#' @param data a `data.frame`, 18#' @param subset see `lm`, 19#' @param weights see `lm`, 20#' @param na.action see `lm`, 21#' @param model one of `"fd"`, `"random"`, or `"pooling"`, 22#' @param index the indexes, see [pdata.frame()], 23#' @param R the number of points for the gaussian quadrature, 24#' @param start a vector of starting values, 25#' @param lower the lower bound for the censored/truncated dependent 26#' variable, 27#' @param upper the upper bound for the censored/truncated dependent 28#' variable, 29#' @param objfun the objective function for the fixed effect model (`model = "fd"`, 30#' irrelevant for other values of the `model` argument ): 31#' one of `"lsq"` for least squares (minimise sum of squares of the residuals) 32#' and `"lad"` for least absolute deviations (minimise sum of absolute values 33#' of the residuals), 34#' @param sample `"cens"` for a censored (tobit-like) sample, 35#' `"trunc"` for a truncated sample, 36#' @param \dots further arguments. 37#' @return For `model = "fd"`, an object of class `c("plm", "panelmodel")`, for 38#' `model = "random"` and `model = "pooling"` an object of class `c("maxLik", "maxim")`. 39#' 40#' @export 41#' @importFrom maxLik maxLik 42#' @author Yves Croissant 43#' @references 44#' 45#' \insertRef{HONO:92}{plm} 46#' 47#' @keywords regression 48#' @examples 49#' ## as these examples take a bit of time, do not run them automatically 50#' \dontrun{ 51#' data("Donors", package = "pder") 52#' library("plm") 53#' pDonors <- pdata.frame(Donors, index = "id") 54#' 55#' # replicate Landry/Lange/List/Price/Rupp (2010), online appendix, table 5a, models A and B 56#' modA <- pldv(donation ~ treatment + prcontr, data = pDonors, 57#' model = "random", method = "bfgs") 58#' summary(modA) 59#' modB <- pldv(donation ~ treatment * prcontr - prcontr, data = pDonors, 60#' model = "random", method = "bfgs") 61#' summary(modB) 62#' } 63#' 64# 65# TODO: check if argument method = "bfgs" is needed in example (and why) 66# -> seems strange as it is no direct argument of pldv 67 68pldv <- function(formula, data, subset, weights, na.action, 69 model = c("fd", "random", "pooling"), index = NULL, 70 R = 20, start = NULL, lower = 0, upper = +Inf, 71 objfun = c("lsq", "lad"), sample = c("cens", "trunc"), ...){ 72 73## Due to the eval() construct with maxLik::maxLik we import maxLik::maxLik 74## and re-export it via NAMESPACE as plm::maxLik with a minimal documentation 75## pointing to the original documentation. 76## This way, we can keep the flexibility of eval() [evalutate in parent frame] 77## and can lessen the dependency burden by placing pkg maxLik in 'Imports' 78## rather than 'Depends' in DESCRIPTION. 79 80 # use the plm interface to compute the model.frame 81 sample <- match.arg(sample) 82 model <- match.arg(model) 83 cl <- match.call() 84 mf <- match.call(expand.dots = FALSE) 85 mf <- cl 86 m <- match(c("formula", "data", "subset", "weights", "na.action", "index"), names(mf), 0) 87 mf <- mf[c(1L, m)] 88 mf$model <- NA 89 mf[[1L]] <- as.name("plm") 90 mf <- eval(mf, parent.frame()) 91 formula <- attr(mf, "formula") 92 93 # extract the relevant arguments for maxLik 94 maxl <- cl 95 m <- match(c("print.level", "ftol", "tol", "reltol", 96 "gradtol", "steptol", "lambdatol", "qrtol", 97 "iterlim", "fixed", "activePar", "method"), names(maxl), 0) 98 maxl <- maxl[c(1L, m)] 99 maxl[[1L]] <- as.name("maxLik") 100 101 # The within model -> Bo Honore (1992) 102 if (model == "fd"){ 103 objfun <- match.arg(objfun) 104 # create a data.frame containing y_t and y_{t-1} 105 y <- as.character(formula[[2L]]) 106 y <- mf[[y]] 107 ly <- c(NA, y[1:(length(y) - 1)]) 108 id <- as.integer(index(mf, "id")) 109 lid <- c(NA, id[1:(nrow(mf) - 1)]) 110 keep <- id == lid 111 keep[1L] <- FALSE 112 Y <- data.frame(y, ly) 113 Y <- Y[keep, ] 114 yt <- Y$y 115 ytm1 <- Y$ly 116 # create the matrix of first differenced covariates 117 X <- model.matrix(mf, model = "fd") 118 start <- rep(.1, ncol(X)) 119 names(start) <- colnames(X) 120 if (sample == "trunc"){ 121 if (objfun == "lad") fm <- function(x) abs(x) 122 if (objfun == "lsq") fm <- function(x) x ^ 2 123 psi <- function(a1, a2, b){ 124 fm( (a2 <= b) * a1 + 125 (b > - a1 & b < a2) * (a2 - a1 - b) + 126 (a1 <= - b) * a2 127 ) 128 } 129 } 130 if (sample == "cens"){ 131 if (objfun == "lad"){ 132 psi <- function(a1, a2, b){ 133 (a1 <= pmax(0, - b) & a2 <= pmax(0, b)) * 0 + 134 (! (a1 <= pmax(0, - b) & a2 <= pmax(0, b)) ) * abs(a2 - a1 - b) 135 } 136 } 137 if (objfun == "lsq"){ 138 psi <- function(a1, a2, b){ 139 (a2 <= b) * (a1 ^ 2 - 2 * a1 * (a2 - b)) + 140 (b > - a1 & b < a2) * (a2 - a1 - b) ^ 2 + 141 (a1 <= - b) * (a2 ^ 2 - 2 * a2 * (b + a1)) 142 } 143 } 144 } 145 BO <- function(param){ 146 bdx <- as.numeric(X %*% param) 147 lnl <- - psi(ytm1, yt, bdx) 148 selobs <- (bdx > - ytm1 & bdx < yt) 149 if (objfun == "lsq" && sample == "cens"){ 150 attr(lnl, "gradient") <- - 151 ( (ytm1 > - bdx & yt > bdx) * (- 2 * (yt - ytm1 - bdx)) + 152 (ytm1 > - bdx & yt < bdx) * ( 2 * ytm1) + 153 (ytm1 < - bdx & yt > bdx) * (- 2 * yt) ) * X 154 attr(lnl, "hessian") <- - crossprod( (ytm1 > - bdx & yt > bdx) * X) 155 } 156 lnl 157 } 158 maxl[c("logLik", "start")] <- list(BO, start) 159 result <- eval(maxl, parent.frame()) 160 if (objfun == "lsq" && sample == "cens"){ 161 bdx <- as.numeric((crossprod(t(X), coef(result)))) 162 V4 <- yt ^ 2 * (bdx <= - ytm1) + ytm1 ^ 2 * (yt <= bdx) + 163 (yt - ytm1 - bdx) ^ 2 * (bdx > - ytm1 & bdx < yt) 164 V4 <- crossprod(X, V4 * X) / length(V4) 165 T4 <- crossprod((bdx > - ytm1 & bdx < yt) * X, X) / length(V4) 166 solve_T4 <- solve(T4) 167 vcov <- solve_T4 %*% V4 %*% solve_T4 168 result$vcov <- V4 169 } 170 if (is.null(result$vcov)) result$vcov <- solve(- result$hessian) 171 resid <- yt - as.numeric(crossprod(t(X), coef(result))) 172 result <- list(coefficients = coef(result), 173 vcov = result$vcov, 174 formula = formula, 175 model = mf, 176 df.residual = nrow(X) - ncol(X), 177 residuals = resid, 178 args = list(model = "fd", effect = "individual"), 179 call = cl) 180 class(result) <- c("plm", "panelmodel") 181 } 182 else{ # model != "fd" => cases model = "random" / "pooling" 183 184 # old pglm stuff for the pooling and the random model, with 185 # update to allow upper and lower bonds 186 X <- model.matrix(mf, rhs = 1, model = "pooling", effect = "individual") 187 188 if (ncol(X) == 0L) stop("empty model") 189 y <- pmodel.response(mf, model = "pooling", effect = "individual") 190 id <- attr(mf, "index")[[1L]] 191 192 # The following is the only instance of statmod::gauss.quad, so check for 193 # the package's availability. (We placed 'statmod' in 'Suggests' rather 194 # than 'Imports' so that it is not an absolutely required dependency.) 195 ## Procedure for pkg check for pkg in 'Suggests' as recommended in 196 ## Wickham, R packages (http://r-pkgs.had.co.nz/description.html). 197 if (!requireNamespace("statmod", quietly = TRUE)) { 198 stop(paste("Function 'gauss.quad' from package 'statmod' needed for this function to work.", 199 "Please install it, e.g., with 'install.packages(\"statmod\")"), 200 call. = FALSE) 201 } 202 # compute the nodes and the weights for the gaussian quadrature 203 rn <- statmod::gauss.quad(R, kind = 'hermite') 204 # compute the starting values 205 ls <- length(start) 206 if (model == "pooling"){ 207 K <- ncol(X) 208 if (! ls %in% c(0, K + 1)) stop("irrelevant length for the start vector") 209 if (ls == 0L){ 210 m <- match(c("formula", "data", "subset", "na.action"), names(cl), 0) 211 lmcl <- cl[c(1,m)] 212 lmcl[[1L]] <- as.name("lm") 213 lmcl <- eval(lmcl, parent.frame()) # eval stats::lm() 214 sig2 <- deviance(lmcl) / df.residual(lmcl) 215 sigma <- sqrt(sig2) 216 start <- c(coef(lmcl), sd.nu = sigma) 217 } 218 } 219 else{ # case model != "pooling" and != "fd" => model ="random" 220 if (ls <= 1L){ 221 startcl <- cl 222 startcl$model <- "pooling" 223 startcl$method <- "bfgs" 224 pglmest <- eval(startcl, parent.frame()) # eval pldv() with updated args 225 thestart <- coef(pglmest) 226 if (ls == 1L){ 227 start <- c(thestart, start) 228 } 229 else{ 230 # case ls = 0 231 resid <- y - as.numeric(tcrossprod(X, t(coef(pglmest)[1:ncol(X)]))) 232 eta <- tapply(resid, id, mean)[as.character(id)] 233 nu <- resid - eta 234 start <- c(thestart[1:ncol(X)], sd.nu = sd(nu), sd.eta = sd(eta)) 235 } 236 } 237 } 238 # call to maxLik with the relevant arguments 239 argschar <- function(args){ 240 paste(as.character(names(args)), as.character(args), 241 sep= "=", collapse= ",") 242 } 243 args <- list(param = "start", 244 y = "y", X = "X", id = "id", model = "model", 245 rn = "rn", lower = lower, upper = upper) 246 thefunc <- paste("function(start) lnl.tobit", "(", argschar(args), ")", sep = "") 247 maxl$logLik <- eval(parse(text = thefunc)) 248 maxl$start <- start 249 result <- eval(maxl, parent.frame()) 250 result[c('call', 'args', 'model')] <- list(cl, args, data) 251 } # end cases model = "random" / "pooling" 252 result 253} 254 255 256lnl.tobit <- function(param, y, X, id, lower = 0, upper = +Inf, model = "pooling", rn = NULL){ 257 compute.gradient <- TRUE 258 compute.hessian <- FALSE 259 mills <- function(x) exp(dnorm(x, log = TRUE) - pnorm(x, log.p = TRUE)) 260 O <- length(y) 261 K <- ncol(X) 262 beta <- param[1L:K] 263 sigma <- param[K + 1L] 264 Xb <- as.numeric(crossprod(t(X), beta)) 265 YLO <- (y == lower) 266 YUT <- (y > lower) & (y < upper) 267 YUP <- y == upper 268 if (model == "random"){ 269 R <- length(rn$nodes) 270 seta <- param[K + 2L] 271 } 272 else seta <- 0 273 274 f <- function(i = NA){ 275 result <- numeric(length = length(y)) 276 if (is.na(i)) z <- 0 else z <- rn$nodes[i] 277 e <- (y - Xb - sqrt(2) * seta * z) / sigma 278 result[YLO] <- pnorm( e[YLO], log.p = TRUE) 279 result[YUT] <- dnorm( e[YUT], log = TRUE) - log(sigma) 280 result[YUP] <- pnorm(- e[YUP], log.p = TRUE) 281 result 282 } 283 284 g <- function(i = NA){ 285 if (is.na(i)) z <- 0 else z <- rn$nodes[i] 286 e <- (y - Xb - sqrt(2) * seta * z) / sigma 287 mz <- mills(e) 288 mmz <- mills(- e) 289 gradi <- matrix(0, nrow = nrow(X), ncol = ncol(X) + 1L) 290 gradi[YLO, 1L:K] <- - mz[YLO] * X[YLO, , drop = FALSE] 291 gradi[YLO, K + 1L] <- - e[YLO] * mz[YLO] 292 gradi[YUT, 1L:K] <- e[YUT] * X[YUT, , drop = FALSE] 293 gradi[YUT, K + 1L] <- - (1 - e[YUT] ^ 2) 294 gradi[YUP, 1L:K] <- mmz[YUP] * X[YUP, , drop = FALSE] 295 gradi[YUP, K + 1L] <- e[YUP] * mmz[YUP] 296 if (! is.na(i)){ 297 gradi <- cbind(gradi, NA) 298 gradi[YLO, K + 2L] <- - mz[YLO] * sqrt(2) * z 299 gradi[YUT, K + 2L] <- e[YUT] * sqrt(2) * z 300 gradi[YUP, K + 2L] < - mmz[YUP] * sqrt(2) * z 301 } 302 gradi / sigma 303 } 304 305 h <- function(i = NA, pwnt = NULL){ 306 if (is.na(i)){ 307 z <- 0 308 seta <- 0 309 pw <- 1 310 } 311 else{ 312 z <- rn$nodes[i] 313 pw <- pwnt[[i]] 314 } 315 e <- (y - Xb - sqrt(2) * seta * z) / sigma 316 mz <- mills(e) 317 mmz <- mills(- e) 318 hbb <- hbs <- hss <- numeric(length = nrow(X)) # pre-allocate 319 hbb[YLO] <- - (e[YLO] + mz[YLO]) * mz[YLO] 320 hbs[YLO] <- mz[YLO] * (1 - (e[YLO] + mz[YLO]) * e[YLO]) 321 hss[YLO] <- e[YLO] * mz[YLO] * (2 - (e[YLO] + mz[YLO]) * e[YLO]) 322 hbb[YUT] <- - 1 323 hbs[YUT] <- - 2 * e[YUT] 324 hss[YUT] <- (1 - 3 * e[YUT] ^ 2) 325 hbb[YUP] <- - (- e[YUP] + mmz[YUP]) * mmz[YUP] 326 hbs[YUP] <- - mmz[YUP] * (1 + (mmz[YUP] - e[YUP]) * e[YUP]) 327 hss[YUP] <- - e[YUP] * mmz[YUP] * (2 + (mmz[YUP] - e[YUP]) * e[YUP]) 328 hbb <- crossprod(hbb * X * pw, X) 329 hbs <- apply(hbs * X * pw, 2, sum) # TODO: can use colSums -> faster 330 hss <- sum(hss * pw) 331 H <- rbind(cbind(hbb, hbs), c(hbs, hss)) 332 if (! is.na(i)){ 333 hba <- hsa <- haa <- numeric(length = nrow(X)) 334 hba[YLO] <- - (e[YLO] + mz[YLO]) * mz[YLO] * sqrt(2) * z 335 hsa[YLO] <- mz[YLO] * sqrt(2) * z * (1 - (e[YLO] + mz[YLO]) * e[YLO]) 336 haa[YLO] <- - (e[YLO] + mz[YLO]) * mz[YLO] * 2 * z ^ 2 337 hba[YUT] <- - sqrt(2) * z 338 hsa[YUT] <- - 2 * sqrt(2) * z * e[YUT] 339 haa[YUT] <- - 2 * z ^ 2 340 hba[YUP] <- - (- e[YUP] + mmz[YUP]) * mmz[YUP] * sqrt(2) * z 341 hsa[YUP] <- - mmz[YUP] * sqrt(2) * z * (1 + (- e[YUP] + mmz[YUP]) * e[YUP]) 342 haa[YUP] <- - (- e[YUP] + mmz[YUP]) * mmz[YUP] * 2 * z ^ 2 343 hba <- apply(hba * X * pw, 2, sum) # TODO: can use colSums -> faster 344 haa <- sum(haa * pw) 345 hsa <- sum(hsa * pw) 346 H <- rbind(cbind(H, c(hba, hsa)), c(hba, hsa, haa)) 347 } 348 H / sigma ^ 2 349 } 350 351 if (model == "pooling"){ 352 lnL <- sum(f(i = NA)) 353 if (compute.gradient) attr(lnL, "gradient") <- g(i = NA) 354 if (compute.hessian) attr(lnL, "hessian") <- h(i = NA) 355 } 356 if (model == "random"){ 357 lnPntr <- lapply(1:R, function(i) f(i = i)) 358 lnPnr <- lapply(lnPntr, function(x){ 359 result <- tapply(x, id, sum) 360 ids <- names(result) 361 result <- as.numeric(result) 362 names(result) <- ids 363 result 364 } 365 ) 366 lnPn <- lapply(1:R, function(i) rn$weights[i] * exp(lnPnr[[i]])) 367 lnPn <- log(Reduce("+", lnPn)) - 0.5 * log(pi) 368 lnL <- sum(lnPn) 369 if (compute.gradient || compute.hessian){ 370 glnPnr <- lapply(1:R, function(i) g(i = i)) 371 pwn <- lapply(1:R, function(i) exp(lnPnr[[i]] - lnPn)) 372 pwnt <- lapply(1:R, function(i) pwn[[i]][as.character(id)]) 373 glnPnr2 <- lapply(1:R, function(i) rn$weights[i] * pwnt[[i]] * glnPnr[[i]]) 374 gradi <- Reduce("+", glnPnr2) / sqrt(pi) 375 attr(lnL, "gradient") <- gradi 376 } 377 if (compute.hessian){ 378 hlnPnr <- lapply(1:R, function(i) h(i = i, pwnt = pwnt)) 379 daub <- lapply(1:R, function(i) apply(glnPnr[[i]], 2, tapply, id, sum) * pwn[[i]] * rn$weights[i]) 380 daub <- Reduce("+", daub) / sqrt(pi) 381 DD1 <- - crossprod(daub) 382 DD2 <- lapply(1:R, function(i) rn$weights[i] * hlnPnr[[i]]) 383 DD2 <- Reduce("+", DD2) / sqrt(pi) 384 DD3 <- lapply(1:R, function(i) rn$weights[i] * crossprod(sqrt(pwn[[i]]) * apply(glnPnr[[i]], 2, tapply, id, sum))) 385 DD3 <- Reduce("+", DD3) / sqrt(pi) 386 H <- (DD1 + DD2 + DD3) 387 attr(lnL, "hessian") <- H 388 } 389 } 390 lnL 391} 392 393 394