R Under development (unstable) (2023-12-13 r85679 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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)
> 
>  # 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)
> 
>  ##### redo with aster.default and predict.aster
> 
>  out2 <- aster(x, root, pred, fam, modmat = out$modmat)
>  sout2 <- summary(out2)
> 
>  foo <- match(sort(unique(site)), site)
>  modmat.pred <- out$modmat[foo, , ]
>  origin.pred <- out$origin[foo, ]
> 
>  pout1 <- predict(out2, modmat = modmat.pred, parm.type = "canon")
> 
>  ##### case 1: model = "unco", obj = "unco", parm = "cano" ####
> 
>  fred <- predict(out2, modmat = modmat.pred, parm.type = "canon",
+      se.fit = TRUE)
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out2$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  sally <- matrix(modmat.pred, ncol = length(out2$coef))
> 
>  all.equal(fred$gradient, sally)
[1] TRUE
> 
>  all.equal(fred$fit, as.numeric(origin.pred) + as.numeric(sally %*% out$coef))
[1] TRUE
> 
>  ##### case 1a: same but with amat
> 
>  node.names <- dimnames(out$modmat)[[2]]
>  site.names <- levels(site)
>  amat <- array(0, c(dim(modmat.pred)[1:2], length(site.names)))
>  for (i in seq(along = site.names))
+      amat[i, grep("h", node.names), i] <- 1
> 
>  alfie <- predict(out2, modmat = modmat.pred, parm.type = "canon",
+      se.fit = TRUE, amat = amat)
> 
>  amatmat <- matrix(amat, ncol = dim(amat)[3])
> 
>  all.equal(alfie$fit, as.numeric(t(amatmat) %*% fred$fit))
[1] TRUE
> 
>  all.equal(alfie$gradient, t(amatmat) %*% fred$gradient)
[1] TRUE
> 
>  all.equal(alfie$se.fit, sqrt(diag(alfie$gradient %*% solve(out2$fisher) %*%
+      t(alfie$gradient))))
[1] TRUE
> 
>  ##### case 2: model = "cond", obj = "cond", parm = "cano" ####
>  ##### no test -- same code as case 1
> 
>  ##### case 3: model = "unco", obj = "cond", parm = "cano" ####
> 
>  out3 <- aster(x, root, pred, fam, modmat = out$modmat, type = "cond")
>  sout3 <- summary(out3)
> 
>  fred <- predict(out3, modmat = modmat.pred, parm.type = "canon",
+      se.fit = TRUE)
> 
>  nind <- dim(modmat.pred)[1]
>  nnode <- dim(modmat.pred)[2]
>  ncoef <- dim(modmat.pred)[3]
> 
>  aster:::setfam(fam.default())
> 
>  beta.hat <- out3$coef
>  theta.hat <- as.numeric(sally %*% beta.hat)
>  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 = double(nind * nnode))$phi
> 
>  all.equal(fred$fit, phi.hat)
[1] TRUE
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out3$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  my.gradient <- 0 * fred$gradient
>  epsilon <- 1e-9
>  for (k in 1:ncoef) {
+      beta.epsilon <- beta.hat
+      beta.epsilon[k] <- beta.hat[k] + epsilon
+      theta.epsilon <- as.numeric(sally %*% beta.epsilon)
+      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 = double(nind * nnode))$phi
+      my.gradient[ , k] <- (phi.epsilon - phi.hat) / epsilon
+  }
> 
>  all.equal(fred$gradient, my.gradient, tolerance = sqrt(epsilon))
[1] TRUE
> 
>  alfie <- predict(out3, modmat = modmat.pred, parm.type = "canon",
+      se.fit = TRUE, amat = amat)
> 
>  all.equal(alfie$fit, as.numeric(t(amatmat) %*% fred$fit))
[1] TRUE
> 
>  all.equal(alfie$gradient, t(amatmat) %*% fred$gradient)
[1] TRUE
> 
>  all.equal(alfie$se.fit, sqrt(diag(alfie$gradient %*% solve(out3$fisher) %*%
+      t(alfie$gradient))))
[1] TRUE
> 
>  ##### case 4: model = "cond", obj = "unco", parm = "cano" ####
> 
>  fred <- predict(out2, modmat = modmat.pred, parm.type = "canon",
+      model.type = "cond", se.fit = TRUE)
> 
>  aster:::setfam(fam.default())
> 
>  beta.hat <- out2$coef
>  phi.hat <- as.numeric(origin.pred) + as.numeric(sally %*% beta.hat)
>  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 = double(nind * nnode))$theta
> 
>  all.equal(fred$fit, theta.hat)
[1] TRUE
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out2$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  my.gradient <- 0 * fred$gradient
>  epsilon <- 1e-9
>  for (k in 1:ncoef) {
+      beta.epsilon <- beta.hat
+      beta.epsilon[k] <- beta.hat[k] + epsilon
+      phi.epsilon <- as.numeric(origin.pred) + as.numeric(sally %*% beta.epsilon)
+      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 = double(nind * nnode))$theta
+      my.gradient[ , k] <- (theta.epsilon - theta.hat) / epsilon
+  }
> 
>  all.equal(fred$gradient, my.gradient, tolerance = sqrt(epsilon))
[1] TRUE
> 
>  alfie <- predict(out2, modmat = modmat.pred, parm.type = "canon",
+      model.type = "cond", se.fit = TRUE, amat = amat)
> 
>  all.equal(alfie$fit, as.numeric(t(amatmat) %*% fred$fit))
[1] TRUE
> 
>  all.equal(alfie$gradient, t(amatmat) %*% fred$gradient)
[1] TRUE
> 
>  all.equal(alfie$se.fit, sqrt(diag(alfie$gradient %*% solve(out2$fisher) %*%
+      t(alfie$gradient))))
[1] TRUE
> 
>  ##### case 5: model = "cond", obj = "cond", parm = "mean" ####
> 
>  root.pred <- matrix(1, nind, nnode)
> 
>  fred <- predict(out3, modmat = modmat.pred, parm.type = "mean",
+      model.type = "cond", root = root.pred, x = root.pred)
> 
>  aster:::setfam(fam.default())
> 
>  beta.hat <- out3$coef
>  theta.hat <- as.numeric(sally %*% beta.hat)
>  xi.hat <- .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.hat),
+      ctau = double(nind * nnode))$ctau
> 
>  all.equal(fred, xi.hat)
[1] TRUE
> 
>  fred <- predict(out3, modmat = modmat.pred, parm.type = "mean",
+      model.type = "cond", root = root.pred, x = root.pred, se.fit = TRUE)
> 
>  all.equal(fred$fit, xi.hat)
[1] TRUE
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out3$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  aster:::setfam(fam.default())
> 
>  my.gradient <- 0 * fred$gradient
>  epsilon <- 1e-9
>  for (k in 1:ncoef) {
+      beta.epsilon <- beta.hat
+      beta.epsilon[k] <- beta.hat[k] + epsilon
+      theta.epsilon <- as.numeric(sally %*% beta.epsilon)
+      xi.epsilon <- .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.epsilon),
+          ctau = double(nind * nnode))$ctau
+      my.gradient[ , k] <- (xi.epsilon - xi.hat) / epsilon
+  }
> 
>  all.equal(fred$gradient, my.gradient, tolerance = sqrt(epsilon))
[1] TRUE
> 
>  ##### case 6: model = "unco", obj = "unco", parm = "mean" ####
> 
>  fred <- predict(out2, modmat = modmat.pred, parm.type = "mean",
+      root = root.pred)
> 
>  beta.hat <- out2$coef
> 
>  beta2tau <- function(beta) {
+ 
+      phi <- origin.pred + matrix(sally %*% 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)
+  }
> 
>  aster:::setfam(fam.default())
> 
>  tau.hat <- beta2tau(beta.hat)
> 
>  all.equal(fred, tau.hat)
[1] TRUE
> 
>  fred <- predict(out2, modmat = modmat.pred, parm.type = "mean",
+      root = root.pred, se.fit = TRUE)
> 
>  all.equal(fred$fit, tau.hat)
[1] TRUE
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out2$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  aster:::setfam(fam.default())
> 
>  my.gradient <- 0 * fred$gradient
>  for (k in 1:length(beta.hat)) {
+      beta.epsilon <- beta.hat
+      beta.epsilon[k] <- beta.hat[k] + epsilon
+      tau.epsilon <- beta2tau(beta.epsilon)
+      my.gradient[ , k] <- (tau.epsilon - tau.hat) / epsilon
+  }
> 
>  all.equal(fred$gradient, my.gradient, tolerance = sqrt(epsilon))
[1] TRUE
> 
>  ##### case 7: model = "cond", obj = "unco", parm = "mean" ####
> 
>  fred <- predict(out2, modmat = modmat.pred, parm.type = "mean",
+      model.type = "cond", root = root.pred, x = root.pred)
> 
>  beta.hat <- out2$coef
> 
>  beta2xi <- function(beta) {
+ 
+      phi <- origin.pred + matrix(sally %*% 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
+ 
+      return(ctau)
+  }
> 
>  aster:::setfam(fam.default())
> 
>  xi.hat <- beta2xi(beta.hat)
> 
>  all.equal(fred, xi.hat)
[1] TRUE
> 
>  fred <- predict(out2, modmat = modmat.pred, parm.type = "mean",
+      model.type = "cond", root = root.pred, x = root.pred, se.fit = TRUE)
> 
>  all.equal(fred$fit, xi.hat)
[1] TRUE
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out2$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  aster:::setfam(fam.default())
> 
>  my.gradient <- 0 * fred$gradient
>  for (k in 1:ncoef) {
+      beta.epsilon <- beta.hat
+      beta.epsilon[k] <- beta.hat[k] + epsilon
+      xi.epsilon <- beta2xi(beta.epsilon)
+      my.gradient[ , k] <- (xi.epsilon - xi.hat) / epsilon
+  }
> 
>  all.equal(fred$gradient, my.gradient, tolerance = sqrt(epsilon))
[1] TRUE
> 
>  ##### case 8: model = "unco", obj = "cond", parm = "mean" ####
> 
>  fred <- predict(out3, modmat = modmat.pred, root = root.pred)
> 
>  beta.hat <- out3$coef
> 
>  beta2tau <- function(beta) {
+ 
+      theta <- matrix(sally %*% 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)
+  }
> 
>  aster:::setfam(fam.default())
> 
>  tau.hat <- beta2tau(beta.hat)
> 
>  all.equal(fred, tau.hat)
[1] TRUE
> 
>  fred <- predict(out3, modmat = modmat.pred, root = root.pred, se.fit = TRUE)
> 
>  all.equal(fred$fit, tau.hat)
[1] TRUE
> 
>  all.equal(fred$se.fit, sqrt(diag(fred$gradient %*% solve(out3$fisher) %*%
+      t(fred$gradient))))
[1] TRUE
> 
>  aster:::setfam(fam.default())
> 
>  my.gradient <- 0 * fred$gradient
>  for (k in 1:ncoef) {
+      beta.epsilon <- beta.hat
+      beta.epsilon[k] <- beta.hat[k] + epsilon
+      tau.epsilon <- beta2tau(beta.epsilon)
+      my.gradient[ , k] <- (tau.epsilon - tau.hat) / epsilon
+  }
> 
>  all.equal(fred$gradient, my.gradient, tolerance = sqrt(epsilon))
[1] TRUE
> 
>  ##### HOORAY !!!!! ##### That's it for aster.predict #####
>  ##### now for aster.predict.formula #####
> 
>  ##### case 1: newdata missing
> 
>  pout2 <- predict(out)
> 
>  newdata <- data.frame(site = factor(LETTERS[1:4]))
>  for (v in vars)
+  newdata[[v]] <- 1
>  newdata$root <- 1
>  newdata$foo <- modmat.pred[ , "l2", "foo"]
> 
>  renewdata <- reshape(newdata, varying = list(vars),
+      direction = "long", timevar = "varb", times = as.factor(vars),
+      v.names = "resp")
> 
>  louise <- predict(out, newdata = renewdata, varvar = varb, idvar = id,
+      root = root, se.fit = TRUE)
> 
>  all.equal(louise$modmat, modmat.pred)
[1] TRUE
> 
>  fred <- predict(out2, modmat = modmat.pred, root = root.pred, se.fit = TRUE)
> 
>  all.equal(louise$fit, fred$fit)
[1] TRUE
> 
>  all.equal(louise$se.fit, fred$se.fit)
[1] TRUE
> 
>  foo <- new.env(parent = emptyenv())
>  bar <- suppressWarnings(try(load("predict.rda", foo), silent = TRUE))
>  if (inherits(bar, "try-error")) {
+      save(sout1, sout2, sout3, pout1, pout2, file = "predict.rda")
+  } else {
+      print(all.equal(sout1, foo$sout1))
+      print(all.equal(sout2, foo$sout2))
+      print(all.equal(sout3, foo$sout3))
+      print(all.equal(pout1, foo$pout1))
+      print(all.equal(pout2, foo$pout2))
+  }
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
> 
>  ##### test for global variables #####
> 
>  saves <- c("out", "renewdata", "out2", "modmat.pred", "root.pred", "louise",
+      "fred")
>  rm(list = setdiff(ls(), saves))
>  ls()
[1] "fred"        "louise"      "modmat.pred" "out"         "out2"       
[6] "renewdata"   "root.pred"  
> 
>  louise.too <- predict(out, newdata = renewdata, varvar = varb, idvar = id,
+      root = root, se.fit = TRUE)
>  identical(louise, louise.too)
[1] TRUE
> 
>  fred.too <- predict(out2, modmat = modmat.pred, root = root.pred,
+      se.fit = TRUE)
>  identical(fred, fred.too)
[1] TRUE
> 
>  ##### test of newcoef #####
> 
>  fake <- out2
>  beta.new <- fake$coefficients + rnorm(length(fake$coefficients)) * 0.1
>  fake$coefficients <- beta.new
>  fred.fake <- predict(fake, modmat = modmat.pred, root = root.pred,
+      se.fit = TRUE)
>  fred.new <- predict(out2, modmat = modmat.pred, root = root.pred,
+      se.fit = TRUE, newcoef = beta.new)
>  identical(fred.fake, fred.new)
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   0.28    0.04    0.32