R/lacm.R

In lacm: Latent Autoregressive Count Models

.pairlik <- function(theta, object, Y = object$Y, X = object$X) {
    if (length(Y) != object$nobs | nrow(X) != object$nobs)
        stop("\nWrong dimensions for Y or X")
    beta <- theta[seq_len(object$p)]
    eta <- X %*% beta + object$offset
    phi <- theta[object$p + 1L]    
    tau2 <- theta[object$p + 2L]
    npairs <- (object$nobs - length(object$kweights)) * length(object$kweights)
    llik <- .C("pairlik",
               as.double(eta),
               as.double(phi),
               as.double(tau2), 
               as.integer(Y),
               as.integer(object$nobs),
               as.integer(length(object$kweights)),
               as.integer(object$latent),
               as.double(object$gh$nodes),
               as.double(object$gh$weights),
               as.integer(object$gh.num),
               output = as.double(rep(0.0, npairs)))$output
    if (all(is.finite(llik))) llik
    else (-sqrt(.Machine$double.xmax))
}

.jacob <- function(theta, object, Y = object$Y, X = object$X) {
    if (length(Y) != object$nobs | nrow(X) != object$nobs)
        stop("\nWrong dimensions for Y or X")
    beta <- theta[seq_len(object$p)]
    eta <- X %*% beta + object$offset
    phi <- theta[object$p + 1L]    
    tau2 <- theta[object$p + 2L]
    npairs <- (object$nobs - length(object$kweights)) * length(object$kweights)
    jac <- .C("jacob",
              as.double(eta),
              as.double(phi),
              as.double(tau2),
              as.integer(Y),
              as.integer(object$nobs),
              as.double(X),
              as.integer(object$p),
              as.integer(length(object$kweights)),
              as.integer(object$latent),
              as.double(object$gh$nodes),
              as.double(object$gh$weights),
              as.integer(object$gh.num),
              output = as.double(rep(0.0, npairs * length(theta))))$output
    if (!all(is.finite(jac)))
        jac <- -sqrt(.Machine$double.xmax)
    matrix(jac, nrow = npairs, ncol = object$p + 2, byrow = TRUE)
}

lacm <- function(formula, data, subset, offset, contrasts = NULL, start.theta = NULL, fixed, d = 1, kernel.type = c("Rectangular", "Trapezoidal"), fit = TRUE, gh.num = 20, reltol.opt = 1e-6, opt.method = c("BFGS", "Nelder-Mead"), maxit.opt = 1000, sandwich.lag = NULL, bread.method = c("Outer-product", "Hessian"), ...) {
    ## lines below inherited from 'glm'
    call <- match.call()
    if (missing(data)) 
        data <- environment(formula)
    mf <- match.call(expand.dots = TRUE)
    m <- match(c("formula", "data", "subset", "offset"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- quote(stats::model.frame)
    mf <- eval(mf, parent.frame())
    mt <- attr(mf, "terms")
    Y <- model.response(mf, "any")
    X <- if (!is.empty.model(mt)) 
             model.matrix(mt, mf, contrasts)
         else matrix(, NROW(Y), 0L)
    offset <- as.vector(model.offset(mf))
    if (!is.null(offset)) {
        if (length(offset) != NROW(Y)) 
            stop(gettextf("number of offsets is %d should equal %d (number of observations)", length(offset), NROW(Y)), domain = NA)
    }
    ## end lines inherited from 'glm'
    if (is.null(offset))
        offset <- rep.int(0, NROW(Y))
    ans <- structure(list(nobs = length(Y), p = NCOL(X), d = round(d), npar = NCOL(X) + 2L, Y = Y, X = X, offset = offset, sandwich.lag = sandwich.lag, fit = fit, gh.num = gh.num, call = match.call(), terms = mt), class = "lacm")
    ans$latent <- TRUE
    ilatent <- c(ans$npar - 1L, ans$npar)
    ## fixed parameters
    if (missing(fixed)) {
        ans$fixed <- rep(NA, ans$npar)       
    }
    else {
        if (length(fixed) != ans$npar)
            stop("fixed has a wrong length")
        ## if tau2 is fixed to zero...
        if (!is.na(fixed[ans$npar])) {
            if (fixed[ans$npar] < sqrt(.Machine$double.eps)) {
                ## then also phi is zero
                fixed[ilatent] <- c(0.0, 0.0)
                ans$latent <- FALSE
            }
        }
        ans$fixed <- fixed
    }
    ans$ifree <- is.na(ans$fixed)
    if (sum(ans$ifree) == 0) ## if all parameter fixed
        ans$fit <- FALSE
    theta <- ans$fixed
    ## kernel weights
    kernel.type <- match.arg(kernel.type)
    dseq <- seq_len(ans$d)
    kw <- switch(kernel.type,
                 "Rectangular" = rep(1.0, ans$d),
                 "Trapezoidal" = c(rep(1.0, ans$d - 1), seq(1.0, 0.0, length = ans$d + 2))[-(2 * d + 1L)])
    ans$kweights <-  kw / sum(kw)  
    ##  starting values
    if (is.null(start.theta)) {
        ## see Davis, Dunsmuir and Wang (1999)
        mod0 <- glm.fit(X, Y, family = poisson())
        mu <- fitted(mod0)
        res <- Y - mu        
        s2 <- sum(res ^ 2 - mu, na.rm = TRUE) / sum(mu ^ 2, na.rm = TRUE)
        s2 <- ifelse(s2 > 0.05, s2, 0.05)
        tau2 <- log(s2 + 1)   
        rho <- sum(res[-1L] * res[-ans$nobs], na.rm = TRUE) / (s2 * sum(mu[-1L] * mu[-ans$nobs], na.rm = TRUE))
        phi <- log(s2 * rho + 1) / tau2
        phi <- ifelse(abs(phi) < 0.95, abs(phi), 0.95) * sign(phi)
        start.theta <- c(coef(mod0)[1L] - 0.5 * tau2, coef(mod0)[-1L], phi, tau2)
    }
    ans$start.theta <- start.theta
    theta[ans$ifree] <- ans$start.theta[ans$ifree]
    names(ans$start.theta) <- c(colnames(X), "phi", "tau2")
    
    ans$objfun <- function(theta.free) {
        theta[ans$ifree] <- theta.free
        if (ans$latent) {
            ## check parameter space
            if (abs(theta[ans$p + 1L]) >= 0.99) return (NA)
            if (theta[ans$p + 2L] < 1e-4) return (NA)
        }
        ## and go
        sum(.pairlik(theta, ans) * ans$kweights)
    }
    ans$grad <- function(theta.free) {
        theta[ans$ifree] <- theta.free
        if (ans$latent) {
            ## check parameter space
            if (abs(theta[ans$p + 1L]) >= 0.99) return (NA)
            if (theta[ans$p + 2L] < 1e-4) return (NA)
        }
        ## and go
        gr <- .colSums(.jacob(theta, ans) * ans$kweights, m =  (ans$nobs - length(ans$kweights)) * length(ans$kweights), n = length(theta))
        gr[ans$ifree]
    }
    if (ans$gh.num <= 0)
        stop("The number of Gauss-Hermite nodes (gh.num) must be a positive integer")
    ## Gauss-Hermite nodes and weights
    ans$gh <-  gauss.quad(ans$gh.num, kind = "hermite")
    if (!ans$fit) {
        ans$plik <- ans$objfun(theta[ans$ifree])
    }
    else {
        ## compute maximum pairwise likelihood estimates
        ans$opt.method <- match.arg(opt.method)
        pl.fit <- try(optim(par = theta[ans$ifree], fn = ans$objfun, gr = ans$grad, control = list(fnscale = -1, reltol = reltol.opt, maxit = maxit.opt), method = ans$opt.method))
        ans$convergence <- pl.fit$convergence
        if (inherits(pl.fit, "try-error"))
            return (pl.fit)
        else {
            theta[ans$ifree] <- pl.fit$par
            ans$plik <- pl.fit$value
        }
    }
    ans$theta <- theta
    names(ans$theta) <- names(ans$start.theta)
    ## variance and CLIC
    ans$jacobian <- .jacob(theta, ans)[ ,ans$ifree]
    bread.method <- match.arg(bread.method)
    ans$outer.product <- identical(bread.method, "Outer-product")
    ans$H <- .bread.lacm(ans, ...)
    if (is.null(ans$sandwich.lag)) 
        ans$sandwich.lag <- floor(10 * log10(ans$nobs))
    ans$J <- .meat.lacm(ans, ...)
    ans$vcov <- ans$H %*% ans$J %*%  ans$H / ans$nobs
    ans$CLIC <- - 2 * ans$plik + 2 * sum(diag(ans$H %*% ans$J))
    colnames(ans$vcov) <- rownames(ans$vcov) <- names(ans$start.theta)[ans$ifree]
    
    class(ans) <- "lacm"    
    ans
}

## this is exactly the same of print.lm
print.lacm <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {
    cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "")
    if (length(coef(x))) {
        cat("Coefficients:\n")
        print.default(format(coef(x), digits = digits), print.gap = 2L, quote = FALSE)
        cat("\n")
    }
    else cat("No coefficients\n")
    cat("\n")
    invisible(x)
}

coef.lacm <- function(object, ...) {
    object$theta
}

.bread.lacm <- function(x, ...) {
    if (x$outer.product) {
        wjacob <- x$jacobian * sqrt(x$kweights)
        allcrossprods <- vapply(1:nrow(wjacob), function(i) tcrossprod(wjacob[i, ]), matrix(0.0, nrow = ncol(wjacob), ncol = ncol(wjacob)))
        H <- rowSums(allcrossprods, dims = 2) / x$nobs
    }
    else{
        H <- jacobian(x$grad, x$theta[x$ifree], method = "simple") / x$nobs
    }
    solve(H)
}

.meat.lacm <- function(x, ...) {
    ## aggregated estimating equations
    id <- rep(1:(x$nobs - length(x$kweights)), each = length(x$kweights))
    efun <- aggregate(x$jacobian  * x$kweights, list(id), sum)[ ,-1]
    efun <- as.matrix(efun)
    colnames(efun) <- names(x$start.theta[x$ifree])
    ## go
    acf.es <- acf(efun, lag.max = x$sandwich.lag, type = "cov", plot = FALSE, demean = FALSE)[[1]]
    vmat <- acf.es[1, , ]
    autocov <- acf.es[-1, , ]
    wb <- (1 - (1:x$sandwich.lag) / x$sandwich.lag)
    if (is.array(autocov))
        wsum <- apply(autocov, c(2,3), function(x) sum(x * wb))
    else
        wsum <- sum(autocov * wb)
    vmat + wsum + t(wsum)
}

vcov.lacm <- function(object, ...) {
    object$vcov
}

CLIC <- function(object, ...) {
    object$CLIC
}

summary.lacm <- function(object, ...) {
    cf <- object$theta
    se <- rep(NA, length(cf))
    se[object$ifree] <- sqrt(diag(object$vcov))
    cf <- cbind(cf, se, cf/se, 2 * pnorm(-abs(cf/se)))
    colnames(cf) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
    rownames(cf) <- names(object$theta)
    object$coefficients <- cf 
    class(object) <- "summary.lacm"
    object
}

print.summary.lacm <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:", deparse(x$call, width.cutoff = floor(getOption("width") * 0.85)), "", sep = "\n")
    cat("Pairwise likelihood order:", x$d, "\n")
    cat("\nCoefficients:\n")
    printCoefmat(x$coefficients, digits = digits, signif.legend = TRUE)
    cat("\nLog pairwise likelihood: ", formatC(x$plik, digits = max(5L, digits + 1L)), ",  CLIC: ", format(x$CLIC, digits = max(4L, digits + 1L)), "\n", sep = "")
    if (x$fit && x$convergence != 0) cat("\nWarning: Model did not converge\n")
    if (!x$fit) cat("\nModel not fitted\n")
    invisible(x)
}

simulate.lacm <- function(object, nsim = 1, seed = NULL, ...) {
    ## lines below inherited from 'simulate.lm'
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
        runif(1)
    if (is.null(seed)) 
        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
    else {
        R.seed <- get(".Random.seed", envir = .GlobalEnv)
        set.seed(seed)
        RNGstate <- structure(seed, kind = as.list(RNGkind()))
        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
    }
    X <- object$X
    beta <- object$theta[seq_len(object$p)]
    phi <- object$theta[object$npar - 1L]
    tau2 <- object$theta[object$npar]

    simone <- function() {
        u <- arima.sim(model = list(ar = phi), n = NROW(X), sd = sqrt(tau2 * (1 - phi ^ 2)))
        rpois(n = NROW(X), lambda = exp(X %*% beta + object$offset + u))
    }
    sims <- replicate(nsim, simone())
    sims <- as.data.frame(sims)
    colnames(sims) <- paste("sim_", seq_len(nsim), sep = "")
    sims
}