swh:1:snp:bbee81fcdc4b36c131a8db323aa6b1ea43209e9a
Tip revision: 1cc05f63437c1a519a5bdc24b3cc669980d81d48 authored by Charles J. Geyer on 13 May 2019, 18:20:03 UTC
version 1.0-3
version 1.0-3
Tip revision: 1cc05f6
gradmat.R
library(aster)
# needed because of the change in R function "sample" in R-devel
suppressWarnings(RNGversion("3.5.2"))
set.seed(42)
nind <- 25
vars <- c("l2", "l3", "f2", "f3", "h2", "h3")
pred <- c(0, 1, 1, 2, 3, 4)
fam <- c(1, 1, 1, 1, 3, 3)
length(pred) == length(fam)
nnode <- length(pred)
theta <- matrix(0, nind, nnode)
root <- matrix(1, nind, nnode)
x <- raster(theta, pred, fam, root)
dimnames(x) <- list(NULL, vars)
data <- as.data.frame(x)
site <- factor(sample(LETTERS[1:4], nind, replace = TRUE))
foo <- rnorm(nind)
data <- data.frame(x, site = site, foo = foo, root = 1)
redata <- reshape(data, varying = list(vars),
direction = "long", timevar = "varb", times = as.factor(vars),
v.names = "resp")
out <- aster(resp ~ foo + site + varb, pred, fam, varb, id, root,
data = redata)
sout1 <- summary(out, show.graph = TRUE)
out2 <- aster(x, root, pred, fam, modmat = out$modmat)
sout2 <- summary(out2)
out3 <- aster(x, root, pred, fam, modmat = out$modmat, type = "cond")
sout3 <- summary(out3)
foo <- new.env(parent = emptyenv())
bar <- suppressWarnings(try(load("gradmat.rda", foo), silent = TRUE))
if (inherits(bar, "try-error")) {
save(sout1, sout2, sout3, file = "gradmat.rda")
} else {
print(all.equal(sout1, foo$sout1))
print(all.equal(sout2, foo$sout2))
print(all.equal(sout3, foo$sout3))
}
foo <- match(sort(unique(site)), site)
modmat.pred <- out$modmat[foo, , ]
##### case 1: eta = theta, zeta = phi
beta.hat <- out3$coef
modmat.pred.mat <- matrix(modmat.pred, ncol = length(beta.hat))
theta.hat <- matrix(modmat.pred.mat %*% beta.hat, nrow = dim(modmat.pred)[1])
nind <- dim(modmat.pred)[1]
nnode <- dim(modmat.pred)[2]
ncoef <- dim(modmat.pred)[3]
aster:::setfam(fam.default())
phi.hat <- .C(aster:::C_aster_theta2phi,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = as.double(theta.hat),
phi = matrix(as.double(0), nind, nnode))$phi
my.phi.hat <- theta.hat
my.phi.hat[ , 4] <- my.phi.hat[ , 4] - log(exp(exp(theta.hat[ , 6])) - 1)
my.phi.hat[ , 3] <- my.phi.hat[ , 3] - log(exp(exp(theta.hat[ , 5])) - 1)
my.phi.hat[ , 2] <- my.phi.hat[ , 2] - log(1 + exp(theta.hat[ , 4]))
my.phi.hat[ , 1] <- my.phi.hat[ , 1] - log(1 + exp(theta.hat[ , 3]))
my.phi.hat[ , 1] <- my.phi.hat[ , 1] - log(1 + exp(theta.hat[ , 2]))
all.equal(my.phi.hat, phi.hat)
gradmat <- 0 * modmat.pred
storage.mode(gradmat) <- "double"
gradmat <- .C(aster:::C_aster_D_beta2theta2phi,
nind = as.integer(nind),
nnode = as.integer(nnode),
ncoef = as.integer(ncoef),
pred = as.integer(pred),
fam = as.integer(fam),
theta = as.double(theta.hat),
modmat = as.double(modmat.pred),
gradmat = gradmat)$gradmat
my.gradmat <- 0 * gradmat
epsilon <- 1e-9
for (k in 1:ncoef) {
beta.epsilon <- beta.hat
beta.epsilon[k] <- beta.hat[k] + epsilon
theta.epsilon <- matrix(modmat.pred.mat %*% beta.epsilon, nrow = nind)
phi.epsilon <- .C(aster:::C_aster_theta2phi,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = as.double(theta.epsilon),
phi = matrix(as.double(0), nind, nnode))$phi
my.gradmat[ , , k] <- (phi.epsilon - phi.hat) / epsilon
}
all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon))
##### case 2: eta = phi, zeta = theta
beta.hat <- out2$coef
phi.hat <- matrix(modmat.pred.mat %*% beta.hat, nrow = nind)
theta.hat <- .C(aster:::C_aster_phi2theta,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
phi = as.double(phi.hat),
theta = matrix(as.double(0), nind, nnode))$theta
gradmat <- .C(aster:::C_aster_D_beta2phi2theta,
nind = as.integer(nind),
nnode = as.integer(nnode),
ncoef = as.integer(ncoef),
pred = as.integer(pred),
fam = as.integer(fam),
theta = as.double(theta.hat),
modmat = as.double(modmat.pred),
gradmat = gradmat)$gradmat
my.gradmat <- 0 * gradmat
epsilon <- 1e-9
for (k in 1:ncoef) {
beta.epsilon <- beta.hat
beta.epsilon[k] <- beta.hat[k] + epsilon
phi.epsilon <- matrix(modmat.pred.mat %*% beta.epsilon, nrow = nind)
theta.epsilon <- .C(aster:::C_aster_phi2theta,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
phi = as.double(phi.epsilon),
theta = matrix(as.double(0), nind, nnode))$theta
my.gradmat[ , , k] <- (theta.epsilon - theta.hat) / epsilon
}
all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon))
##### case 3: eta = phi, zeta = tau
root.pred <- matrix(1, nind, nnode)
beta.hat <- out2$coef
beta2tau <- function(beta) {
phi <- matrix(modmat.pred.mat %*% beta, nrow = nind)
theta <- .C(aster:::C_aster_phi2theta,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
phi = as.double(phi),
theta = matrix(as.double(0), nind, nnode))$theta
ctau <- .C(aster:::C_aster_theta2ctau,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = as.double(theta),
ctau = double(nind * nnode))$ctau
tau <- .C(aster:::C_aster_ctau2tau,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
root = as.double(root.pred),
ctau = as.double(ctau),
tau = double(nind * nnode))$tau
return(tau)
}
tau.hat <- beta2tau(beta.hat)
gradmat <- .C(aster:::C_aster_D_beta2phi2tau,
nind = as.integer(nind),
nnode = as.integer(nnode),
ncoef = as.integer(ncoef),
pred = as.integer(pred),
fam = as.integer(fam),
beta = as.double(beta.hat),
root = as.double(root.pred),
origin = rep(as.double(0), nind * nnode),
modmat = as.double(modmat.pred),
gradmat = gradmat)$gradmat
my.gradmat <- 0 * gradmat
epsilon <- 1e-9
for (k in 1:ncoef) {
beta.epsilon <- beta.hat
beta.epsilon[k] <- beta.hat[k] + epsilon
tau.epsilon <- beta2tau(beta.epsilon)
my.gradmat[ , , k] <- (tau.epsilon - tau.hat) / epsilon
}
all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon))
##### case 4: eta = theta, zeta = tau
beta.hat <- out3$coef
beta2tau <- function(beta) {
theta <- matrix(modmat.pred.mat %*% beta, nrow = nind)
ctau <- .C(aster:::C_aster_theta2ctau,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = as.double(theta),
ctau = double(nind * nnode))$ctau
tau <- .C(aster:::C_aster_ctau2tau,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
root = as.double(root.pred),
ctau = as.double(ctau),
tau = double(nind * nnode))$tau
return(tau)
}
tau.hat <- beta2tau(beta.hat)
gradmat <- .C(aster:::C_aster_D_beta2theta2tau,
nind = as.integer(nind),
nnode = as.integer(nnode),
ncoef = as.integer(ncoef),
pred = as.integer(pred),
fam = as.integer(fam),
beta = as.double(beta.hat),
root = as.double(root.pred),
modmat = as.double(modmat.pred),
gradmat = gradmat)$gradmat
my.gradmat <- 0 * gradmat
for (k in 1:ncoef) {
beta.epsilon <- beta.hat
beta.epsilon[k] <- beta.hat[k] + epsilon
tau.epsilon <- beta2tau(beta.epsilon)
my.gradmat[ , , k] <- (tau.epsilon - tau.hat) / epsilon
}
all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon))