R/estimateCt.R

In roualdes/spiders: Fits Predator Preferences Model

##' @title estimate hypothesis c_st = c_t
##'
##' @description estimates parameters from hypothesis lambda_t = c_t * gamma_t
##'
##' @details There are S*T + T free parameters under this hypothesis.
##' 
##' @param Xdst matrix of sums of number of eaten prey species s during occurrence t; rows indexed by time, and cols indexed by prey species, TxS
##' @param Ydst matrix sum of number of caught prey species s during occurrence t; rows indexed by time, and cols indexed by prey species, TxS
##' @param J vector of predators caught in each time period
##' @param I vector of number of days all traps were left out in a given time period
##' @param EM boolean; whether or not EM algorithm is used
##' @param em_maxiter integer specifying max number of EM iterations
##' @param BALANCED boolean; whether or not data are BALANCED
estCt <- function(Xdst, Ydst, J, I, EM, em_maxiter, BALANCED) {

    ## some numbers
    S <- ncol(Xdst); s <- seq_len(S)
    T <- nrow(Xdst); t <- seq_len(T)
    ST <- S*T

    if (EM) {
        ## avoiding estimated zeros
        pen <- matrix(0, T, S)
        zos <- which(Xdst==0, arr.ind=TRUE)
        if (length(zos)>0) {
            pen[zos] <- 1
            pen[zos] <- pen[zos]/(J[zos[,1]]+1) # add imaginary observation
            pen <- 0.001*pen                    # weight it
        }
        
        ## initialize some values
        em_iter <- 1
        ## cHat <- cHat_old <- rep(0.5, length(J)) # runif(length(J))
        init <- est1(Xdst, Ydst, J, I, EM, em_maxiter, BALANCED)
        gammaHat <- gammaHat_old <- init$gamma
        cHat <- cHat_old <- as.vector(sumSp(Xdst) / (J*sumSp(gammaHat)) + sumSp(pen))

        ## iterate EM
        while ( TRUE ) {
            
            ## expected value of Xjst
            lambda <- cHat*gammaHat
            elambda <- exp(-lambda)
            EX <- lambda/(1-elambda)

            ## convenience
            ZEX <- Xdst*EX + pen        # pull ZEX away from zero

            ## iterate only once simultaneous eqs
            gammaHat <- (ZEX + Ydst) / (cHat*J + I)
            cHat <- sumSp(ZEX) / (J*sumSp(gammaHat))
            
            ## check convergence of EM
            if ( converged(cHat, cHat_old) &&
                converged(gammaHat, gammaHat_old) ) {
                break
            }
            
            ## if not converged, store updated estimates
            cHat_old <- cHat
            gammaHat_old <- gammaHat
            ## print(sprintf('EM iteration %d found values: c = %f', em_iter, cHat))
            em_iter <- em_iter+1
            
            ## limit iterations
            if ( em_iter > em_maxiter ) {
                stop(sprintf("estCt: max EM iterations, %d, reached. Please adjust accordingly.", em_maxiter))
            }
            
        }

        ## calc standard error with est params
        ## SE <- seEM(NULL, gammaHat, cHat, Xdst, Ydst, J, I)
        Info <- diag(ST+T)             # initialize information matrix
        g2 <- gammaHat^2               # gamma^2
        l <- cHat*gammaHat
        expl <- exp(l); tmp <- J*expl/(expl-1)^2 # a common term

        ## fill Info with second derivatives
        diag(Info)[-t] <- unlist(Ydst/g2 + cHat^2*tmp) # gamma
        diag(Info)[t] <- sumSp(tmp*g2)                 # c
        for ( i in t ) {           # fill in off diags; upper tri only
            Info[i,-t][seq(0, (S-1)*T, by=T)+i] <- unlist(-J[i]/(expl[i,] - 1) + l[i,]*tmp[i,])
        }
        lowmat <- lower.tri(Info)
        Info[lowmat] <- t(Info)[lowmat]       # make symmetric from upper tri
        tryCatch(var <- solve(Info),
                 error=function(e) {
                     print("Variances not calculated; system is singular.")
                     var <- NULL
                 })
        
        ## calc log-lik with est params
        loglik <- llEM(Xdst, Ydst, NA, gammaHat, J, I, cHat)
        
        list(c=cHat, gamma=as.matrix(gammaHat), em_iters=em_iter,
             ll=loglik, var=var)
    } else {
        
        ## some numbers
        XYdst <- Xdst + Ydst
        stXdst <- sumST(Xdst)
        iter <- 1; maxiter <- 500

        ## not sure this is the right spot for these checks
        ## ensure J & I have dimension T or 1
        if ( length(J) != T ) {
            stop("J indexed oddly says est0.")
        }
        if ( length(I) != T ) {
            stop("I indexed oddly says est0.")
        }

        ## initialize some values
        ## cHat <- cHat_old <- rep(0.5, length(J)) # runif(length(J))
        ## gammaHat <- gammaHat_old <- XYdst / (J*cHat + I)
        gammaHat <- gammaHat_old <- XYdst / (J + I)

        ## iteratively update; relies on concavity of log-lik
        while ( TRUE ) {

            ## update parameters
            cHat <- sumSp(Xdst) / (J*sumSp(gammaHat))
            gammaHat <- XYdst / (J*cHat + I) # row-wise division

            ## check convergence
            if ( converged(gammaHat, gammaHat_old) &&
                converged(cHat, cHat_old) ) {
                break
            }

            ## if not converged, update estimates for next iteration
            gammaHat_old <- gammaHat
            cHat_old <- cHat
            iter <- iter+1

            ## limit iterations
            if ( iter > maxiter ) {
                stop(sprintf("estCt: %d not sufficient iterations for simultaneous equations.", maxiter))
            }
            
        }

        ## calc standard error with est params
        ## SE <- se(NULL, gammaHat, cHat, Xdst, Ydst, J, I)
        Info <- diag(ST+T)              # initialize information matrix

        ## fill Info with second derivatives
        diag(Info)[-t] <- unlist(XYdst/gammaHat^2) # gamma
        diag(Info)[t] <- sumSp(Xdst)/cHat^2        # c
        for ( i in t ) {              # fill off diags; upper tri only
            Info[i,-t][seq(0, (S-1)*T, by=T)+i] <- J[i]
        }
        lowmat <- lower.tri(Info)
        Info[lowmat] <- t(Info)[lowmat]       # make symmetric from upper tri

        ## calc log-lik with est params
        loglik <- ll(Xdst, Ydst, NA, gammaHat, J, I, cHat)
        
        list(gamma=as.matrix(gammaHat), c=cHat, iters=iter,
             ll=loglik, var=solve(Info))
    }
}