inst/doc/delta.R

In aster: Aster Models

### R code from vignette source 'delta.Rnw'

###################################################
### code chunk number 1: options-width
###################################################
options(keep.source = TRUE, width = 80)


###################################################
### code chunk number 2: libraries
###################################################
library(aster)
library(numDeriv)


###################################################
### code chunk number 3: data
###################################################
data(sim)


###################################################
### code chunk number 4: fit.aster
###################################################
aout <- aster(resp ~ varb + 0 + z1 + z2 + I(z1^2) + I(z1 * z2) + I(z2^2),
    pred, fam, varb, id, root, data = redata)
summary(aout)


###################################################
### code chunk number 5: landscape.maximum.function
###################################################
foo <- function(beta) {
    A <- matrix(NaN, 2, 2)
    A[1, 1] <- beta["I(z1^2)"]
    A[2, 2] <- beta["I(z2^2)"]
    A[1, 2] <- beta["I(z1 * z2)"] / 2
    A[2, 1] <- beta["I(z1 * z2)"] / 2
    b <- rep(NaN, 2)
    b[1] <- beta["z1"]
    b[2] <- beta["z2"]
    solve(-2 * A, b)
}


###################################################
### code chunk number 6: try.one
###################################################
foo(aout$coef)


###################################################
### code chunk number 7: new.way
###################################################
Sigma <- vcov(aout)
B <- jacobian(foo, aout$coef)
V <- B %*% Sigma %*% t(B)
V


###################################################
### code chunk number 8: debug-vcov
###################################################
all.equal(Sigma, solve(aout$fisher), check.attributes = FALSE)


###################################################
### code chunk number 9: jacobian.check
###################################################
beta <- aout$coef
A <- matrix(NaN, 2, 2)
A[1, 1] <- beta["I(z1^2)"]
A[2, 2] <- beta["I(z2^2)"]
A[1, 2] <- beta["I(z1 * z2)"] / 2
A[2, 1] <- beta["I(z1 * z2)"] / 2
b <- rep(NaN, 2)
b[1] <- beta["z1"]
b[2] <- beta["z2"]
jack <- matrix(0, nrow = nrow(B), ncol = ncol(B))
# d b / d beta["z1"]
i <- names(beta) == "z1"
jack[ , i] <- solve(-2 * A, c(1, 0))
# d b / d beta["z2"]
i <- names(beta) == "z2"
jack[ , i] <- solve(-2 * A, c(0, 1))
# d A / d beta["I(z1^2)"]
dA <- matrix(0, 2, 2)
dA[1, 1] <- 1
i <- names(beta) == "I(z1^2)"
jack[ , i] <- (1 / 2) * solve(A) %*% dA %*% solve(A) %*% b
# d A / d beta["I(z2^2)"]
dA <- matrix(0, 2, 2)
dA[2, 2] <- 1
i <- names(beta) == "I(z2^2)"
jack[ , i] <- (1 / 2) * solve(A) %*% dA %*% solve(A) %*% b
# d A / d beta["I(z1 * z2)"]
dA <- matrix(0, 2, 2)
dA[1, 2] <- 1 / 2
dA[2, 1] <- 1 / 2
i <- names(beta) == "I(z1 * z2)"
jack[ , i] <- (1 / 2) * solve(A) %*% dA %*% solve(A) %*% b
all.equal(jack, B)


###################################################
### code chunk number 10: clear
###################################################
rm(list = ls())


###################################################
### code chunk number 11: reaster
###################################################
lout <- suppressWarnings(try(load("rout.rda"), silent = TRUE))
if (inherits(lout, "try-error")) {
    data(grey_cloud_2015)
    modmat.sire <- model.matrix(~ 0 + fit:paternalID, redata)
    modmat.dam <- model.matrix(~ 0 + fit:maternalID, redata)
    modmat.siredam <- cbind(modmat.sire, modmat.dam)
    rout.time <- system.time(
        rout <- reaster(resp ~ fit + varb,
            list(parental = ~ 0 + modmat.siredam, block = ~ 0 + fit:block),
            pred, fam, varb, id, root, data = redata)
    )
    save(rout, rout.time, file = "rout.rda")
}


###################################################
### code chunk number 12: rout.time
###################################################
rout.time.seconds <- rout.time["elapsed"] %% 60
rout.time.minutes <- rout.time["elapsed"] %/% 60


###################################################
### code chunk number 13: vcov
###################################################
Sigma <- vcov(rout, standard.deviation = FALSE, re.too = TRUE)


###################################################
### code chunk number 14: map-factory
###################################################
map.factory <- function(rout) {
    stopifnot(inherits(rout, "reaster"))
    aout <- rout$obj
    stopifnot(inherits(aout, "aster"))
    nnode <- ncol(aout$x)
    nind <- nrow(aout$x)
    fixed <- rout$fixed
    random <- rout$random
    if (nnode == 4) {
        is.subsamp <- rep(FALSE, 4)
    } else if (nnode == 5) {
        is.subsamp <- c(FALSE, FALSE, FALSE, TRUE, FALSE)
    } else stop("can only deal with graphs for individuals with 4 or 5 nodes",
        "\nand graph is linear, and subsampling arrow is 4th of 5")
    # fake object of class aster
    randlab <- unlist(lapply(rout$random, colnames))
    include.random <- grepl("paternalID", randlab, fixed = TRUE)
    fake.out <- aout
    fake.beta <- with(rout, c(alpha, b[include.random]))
    modmat.random <- Reduce(cbind, random)
    stopifnot(ncol(modmat.random) == length(rout$b))
    # never forget drop = FALSE in programming R
    modmat.random <- modmat.random[ , include.random, drop = FALSE]
    fake.modmat <- cbind(fixed, modmat.random)
    # now have to deal with objects of class aster (as opposed to reaster)
    # thinking model matrices are three-way arrays.
    stopifnot(prod(dim(aout$modmat)[1:2]) == nrow(fake.modmat))
    fake.modmat <- array(as.vector(fake.modmat),
        dim = c(dim(aout$modmat)[1:2], ncol(fake.modmat)))
    fake.out$modmat <- fake.modmat
    nparm <- length(rout$alpha) + length(rout$b) + length(rout$nu)
    is.alpha <- 1:nparm %in% seq_along(rout$alpha)
    is.bee <- 1:nparm %in% (length(rout$alpha) + seq_along(rout$b))
    is.nu <- (! (is.alpha | is.bee))
    # figure out individuals from each family
    m <- rout$random$parental
    dads <- grep("paternal", colnames(m))
    # get family, that is, paternalID or grandpaternalID as the case may be
    fams <- colnames(m)[dads] |> sub("^.*ID", "", x = _)
    # drop maternal effects columns (if any)
    m.dads <- m[ , dads, drop = FALSE]
    # make into 3-dimensional array, like obj$modmat
    m.dads <- array(m.dads, c(nind, nnode, ncol(m.dads)))
    # only keep fitness node
    # only works for linear graph
    m.dads <- m.dads[ , nnode, ]
    # redefine dads as families of individuals
    stopifnot(as.vector(m.dads) %in% c(0, 1))
    stopifnot(rowSums(m.dads) == 1)
    # tricky, only works because each row of m.dads
    # is indicator vector of family,
    # so we are multiplying family number by zero or one
    dads <- drop(m.dads %*% as.integer(fams))
    # find one individual in each family
    sudads <- sort(unique(dads))
    which.ind <- match(sudads, dads)
    function(alphabeenu) {
        stopifnot(is.numeric(alphabeenu))
        stopifnot(is.finite(alphabeenu))
        stopifnot(length(alphabeenu) == nparm)
        alpha <- alphabeenu[is.alpha]
        bee <- alphabeenu[is.bee]
        nu <- alphabeenu[is.nu]
        fake.beta <- c(alpha, bee[include.random])
        fake.out$coefficients <- fake.beta
        pout <- predict(fake.out, model.type = "conditional",
            is.always.parameter = TRUE)
        xi <- matrix(pout, ncol = nnode)
        xi <- xi[ , ! is.subsamp, drop = FALSE]
        mu <- apply(xi, 1, prod)
        mu <- mu[which.ind]
        names(mu) <- paste0("PID",
            formatC(sudads, format="d", width=3, flag="0"))
        return(mu)
    }
}


###################################################
### code chunk number 15: map.factory.try
###################################################
alphabeenu <- with(rout, c(alpha, b, nu))
map <- map.factory(rout)
mu.hat <- map(alphabeenu)
mu.hat


###################################################
### code chunk number 16: mu.vcov
###################################################
jack <- jacobian(map, alphabeenu)
Sigma.mu.hat <- jack %*% Sigma %*% t(jack)


###################################################
### code chunk number 17: mf
###################################################
fitness_change <- function(mu) mean(mu * (mu / mean(mu) - 1))
delta.fitness <- fitness_change(mu.hat)
delta.fitness


###################################################
### code chunk number 18: delta.fit.vcov
###################################################
jack <- jacobian(fitness_change, mu.hat)
Sigma.delta.fit <- jack %*% Sigma.mu.hat %*% t(jack)
Sigma.delta.fit


###################################################
### code chunk number 19: delta.fit.se
###################################################
sqrt(drop(Sigma.delta.fit))