https://github.com/cran/aster
Tip revision: aa47935123bfca8a22cbc8345d658d0c1713a289 authored by Charles J. Geyer on 14 December 2023, 15:20:02 UTC
version 1.1-3
version 1.1-3
Tip revision: aa47935
gradmat.Rout.save
R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
>
> library(aster)
Loading required package: trust
>
> # 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)
[1] TRUE
> 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))
+ }
[1] TRUE
[1] TRUE
[1] TRUE
>
> 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)
[1] TRUE
>
> 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))
[1] TRUE
>
> ##### 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))
[1] TRUE
>
> ##### 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))
[1] TRUE
>
> ##### 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))
[1] TRUE
>
>
> proc.time()
user system elapsed
0.225 0.027 0.246