### 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))