R Under development (unstable) (2025-05-30 r88253 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 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)
> 
>  data(radish)
> 
>  pred <- c(0,1,2)
>  fam <- c(1,3,2)
> 
>  ### need object of type aster to find idrop below
> 
>  aout <- aster(resp ~ varb + fit : (Site * Region),
+      pred, fam, varb, id, root, data = radish)
> 
>  ### model matrices for fixed and random effects
> 
>  modmat.fix <- model.matrix(resp ~ varb + fit : (Site * Region), data = radish)
>  modmat.blk <- model.matrix(resp ~ 0 + fit : Block, data = radish)
>  modmat.pop <- model.matrix(resp ~ 0 + fit : Pop, data = radish)
> 
>  rownames(modmat.fix) <- NULL
>  rownames(modmat.blk) <- NULL
>  rownames(modmat.pop) <- NULL
> 
>  idrop <- match(aout$dropped, colnames(modmat.fix))
>  idrop <- idrop[! is.na(idrop)]
>  modmat.fix <- modmat.fix[ , - idrop, drop = FALSE]
> 
>  nfix <- ncol(modmat.fix)
>  nblk <- ncol(modmat.blk)
>  npop <- ncol(modmat.pop)
> 
>  ### reaster
> 
>  rout <- reaster(modmat.fix, list(block = modmat.blk, pop = modmat.pop),
+      pred, fam, radish$varb, radish$id, radish$root, response = radish$resp)
> ## IGNORE_RDIFF_BEGIN
>  summary(rout, standard.deviation = FALSE)

Fixed Effects:
                            Estimate Std. Error  z value Pr(>|z|)    
(Intercept)               -467.24226    1.75183 -266.717   <2e-16 ***
varbFlowers                474.13817    1.75417  270.293   <2e-16 ***
varbFruits                 466.11033    1.76038  264.779   <2e-16 ***
fit:SitePoint Reyes         -0.03620    0.20781   -0.174    0.862    
fit:RegionS                 -0.12249    0.07892   -1.552    0.121    
fit:SiteRiverside:RegionS    0.49930    0.01211   41.223   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Variance Components (P-values are one-tailed):
      Estimate Std. Error z value Pr(>|z|)/2  
block 0.107714   0.048296   2.230     0.0129 *
pop   0.009252   0.005756   1.607     0.0540 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> ## IGNORE_RDIFF_END
> 
>  ### can we figure out idrop from here ???
>  identical(colnames(modmat.fix), dimnames(rout$obj$modmat)[[3]])
[1] TRUE
> 
>  ### make zwz
> 
>  zwz <- makezwz(rout$sigma, c(rout$alpha, rout$c), rout$fixed, rout$random,
+      rout$obj)
> 
>  # objective functions
> 
>  objfun_beenu <- objfun.factory(modmat.fix,
+      list(block = modmat.blk, pop = modmat.pop),
+      radish$resp, pred, fam, radish$root, zwz = zwz,
+      standard = FALSE, deriv = 2)
> 
>  objfun_ceesigma <- objfun.factory(modmat.fix,
+      list(block = modmat.blk, pop = modmat.pop),
+      radish$resp, pred, fam, radish$root, zwz = zwz,
+      standard = TRUE, deriv = 2)
> 
>  alphabeenu <- as.vector(with(rout, c(alpha, b, nu)))
>  alphaceesigma <- as.vector(with(rout, c(alpha, c, sigma)))
> 
>  oout_beenu <- objfun_beenu(alphabeenu)
>  oout_ceesigma <- objfun_ceesigma(alphaceesigma)
> 
>  # check the former
> 
>  is.alpha <- seq_along(alphabeenu) <= length(rout$alpha)
>  is.nu_or_sigma <- seq_along(alphabeenu) > length(rout$alpha) + length(rout$b)
>  is.bee_or_cee <- (! (is.alpha | is.nu_or_sigma))
> 
>  alphabee <- with(rout, c(alpha, b))
>  modmat <- with(rout, Reduce(cbind, random, init = fixed))
>  moo <- mlogl(alphabee, pred, fam, radish$resp, radish$root, modmat, deriv = 2)
>  alpha <- alphabeenu[is.alpha]
>  bee <- alphabeenu[is.bee_or_cee]
>  nu <- alphabeenu[is.nu_or_sigma]
>  zwz <- rout$zwz
> 
>  nrand <- sapply(rout$random, ncol)
>  dee_idx <- rep(seq_along(nrand), times = nrand)
>  dee <- nu[dee_idx]
> 
>  # value of objective function, equation (4.5) in Aster Theory book
>  all.equal(oout_beenu$value, as.numeric(moo$value + sum(bee^2 / dee) / 2.0 +
+      determinant(zwz %*% diag(dee) + diag(length(dee)))$modulus / 2.0))
[1] TRUE
>  all.equal(oout_ceesigma$value, oout_beenu$value)
[1] TRUE
> 
>  # gradient with respect to (alpha, b, nu)
>  # equation (4.12a) in Aster Theory book
>  all.equal(oout_beenu$gradient[is.alpha], moo$gradient[is.alpha])
[1] TRUE
>  # equation (4.12b) in Aster Theory book
>  all.equal(oout_beenu$gradient[is.bee_or_cee], moo$gradient[is.bee_or_cee]
+      + bee / dee)
[1] TRUE
>  # equation (4.12c) in Aster Theory book
>  foompter <- rep(0, length(nu))
>  for (k in seq_along(foompter)) {
+      eek <- as.numeric(dee_idx == k)
+      foompter[k] <- sum(- bee^2 / dee^2 * eek) / 2.0 +
+          sum(diag(solve(zwz %*% diag(dee) + diag(length(dee))) %*%
+              zwz %*% diag(eek))) / 2.0
+  }
>  all.equal(oout_beenu$gradient[is.nu_or_sigma], foompter)
[1] TRUE
>  # value and gradient w. r. t. (alpha, b, nu) now checked
> 
>  # check symmetry of Hessians
>  isSymmetric(oout_beenu$hessian)
[1] TRUE
>  isSymmetric(oout_ceesigma$hessian)
[1] TRUE
> 
>  # on to Hessian w. r. t. (alpha, bee, nu)
>  # equation (4.18a) in Aster Theory book
>  is.alpha.short <- seq_along(alphabee) <= length(alpha)
>  is.bee.short <- (! is.alpha.short)
>  all.equal(oout_beenu$hessian[is.alpha, is.alpha],
+      moo$hessian[is.alpha.short, is.alpha.short])
[1] TRUE
>  # equation (4.18b) in Aster Theory book
>  all.equal(oout_beenu$hessian[is.alpha, is.bee_or_cee],
+      moo$hessian[is.alpha.short, is.bee.short])
[1] TRUE
>  all.equal(oout_beenu$hessian[is.bee_or_cee, is.alpha],
+      moo$hessian[is.bee.short, is.alpha.short])
[1] TRUE
>  # equation (4.18c) in Aster Theory book
>  all.equal(oout_beenu$hessian[is.bee_or_cee, is.bee_or_cee],
+      moo$hessian[is.bee.short, is.bee.short] + diag(1 / dee))
[1] TRUE
>  # equation (4.18d) in Aster Theory book
>  foompter <- matrix(0, length(alpha), length(nu))
>  all.equal(oout_beenu$hessian[is.alpha, is.nu_or_sigma], foompter)
[1] TRUE
>  all.equal(oout_beenu$hessian[is.nu_or_sigma, is.alpha], t(foompter))
[1] TRUE
>  # equation (4.18e) in Aster Theory book
>  foompter <- matrix(0, length(bee), length(nu))
>  for (k in seq_along(nu)) {
+      eek <- as.numeric(dee_idx == k)
+      foompter[ , k] <- (- eek * bee / nu[k]^2)
+  }
>  all.equal(oout_beenu$hessian[is.bee_or_cee, is.nu_or_sigma], foompter)
[1] TRUE
>  all.equal(oout_beenu$hessian[is.nu_or_sigma, is.bee_or_cee], t(foompter))
[1] TRUE
>  # equation (4.18f) in Aster Theory book
>  inv <- solve(zwz %*% diag(dee) + diag(length(dee)))
>  foompter <- matrix(0, length(nu), length(nu))
>  for (k1 in seq_along(nu)) {
+      eek1 <- as.numeric(dee_idx == k1)
+      foompter[k1, k1] <- sum(bee^2 / dee^3 * eek1^2)
+      for (k2 in seq_along(nu)) {
+          eek2 <- as.numeric(dee_idx == k2)
+          foompter[k1, k2] <- foompter[k1, k2] -
+              sum(diag(inv %*% zwz %*% diag(eek1) %*%
+              inv %*% zwz %*% diag(eek2))) / 2.0
+      }
+  }
>  all.equal(oout_beenu$hessian[is.nu_or_sigma, is.nu_or_sigma], foompter)
[1] TRUE
>  # value and gradient and Hessian w. r. t. (alpha, b, nu) now checked
> 
>  # gradient w. r. t. (alpha, c, sigma)
>  cee <- alphaceesigma[is.bee_or_cee]
>  sigma <- alphaceesigma[is.nu_or_sigma]
>  aaa <- sqrt(dee)
>  # equation (4.33a) in Aster Theory book
>  identical(oout_beenu$gradient[is.alpha], moo$gradient[is.alpha.short])
[1] TRUE
>  # equation (4.33b) in Aster Theory book
>  all.equal(oout_ceesigma$gradient[is.bee_or_cee],
+      aaa * moo$gradient[is.bee.short] + cee)
[1] TRUE
>  # same, but from chain rule
>  all.equal(oout_ceesigma$gradient[is.bee_or_cee],
+      oout_beenu$gradient[is.bee_or_cee] * aaa)
[1] TRUE
>  # equation (4.33c) in Aster Theory book
>  foompter <- rep(0, length(sigma))
>  for (k in seq_along(foompter)) {
+      eek <- as.numeric(dee_idx == k)
+      foompter[k] <- sum(oout_beenu$gradient[is.bee_or_cee] * cee * eek) +
+          oout_beenu$gradient[is.nu_or_sigma][k] * 2 * sigma[k]
+  }
>  all.equal(oout_ceesigma$gradient[is.nu_or_sigma], foompter)
[1] TRUE
>  # value and gradient w. r. t. (alpha, c, sigma) now checked
> 
>  # Hessian w. r. t. (alpha, c, sigma)
>  # equation (4.35a) in Aster Theory book
>  identical(oout_beenu$hessian[is.alpha, is.alpha],
+      oout_ceesigma$hessian[is.alpha, is.alpha])
[1] TRUE
>  # equation (4.35b) in Aster Theory book
>  all.equal(oout_ceesigma$hessian[is.alpha, is.bee_or_cee],
+      oout_beenu$hessian[is.alpha, is.bee_or_cee] %*% diag(aaa))
[1] TRUE
>  all.equal(oout_ceesigma$hessian[is.bee_or_cee, is.alpha],
+      diag(aaa) %*% oout_beenu$hessian[is.bee_or_cee, is.alpha])
[1] TRUE
>  # equation (4.35c) in Aster Theory book
>  foompter <- matrix(0, length(alpha), length(sigma))
>  for (k in seq_along(sigma)) {
+      eek <- as.numeric(dee_idx == k)
+      foompter[ , k] <- drop(oout_beenu$hessian[is.alpha, is.bee_or_cee] %*%
+          cbind(eek * cee))
+  }
>  all.equal(oout_ceesigma$hessian[is.alpha, is.nu_or_sigma], foompter)
[1] TRUE
>  all.equal(oout_ceesigma$hessian[is.nu_or_sigma, is.alpha], t(foompter))
[1] TRUE
>  # equation (4.35d) in Aster Theory book
>  all.equal(oout_ceesigma$hessian[is.bee_or_cee, is.bee_or_cee],
+      diag(aaa) %*% moo$hessian[is.bee.short, is.bee.short] %*% diag(aaa)
+      + diag(sum(is.bee.short)))
[1] TRUE
>  # equation (4.35e) in Aster Theory book
>  foompter <- matrix(0, length(bee), length(sigma))
>  for (k in seq_along(sigma)) {
+      eek <- as.numeric(dee_idx == k)
+      moobb <- moo$hessian[is.bee.short, is.bee.short]
+      foompter[ , k] <- drop(diag(aaa) %*% moobb %*% cbind(eek * cee)) +
+          moo$gradient[is.bee.short] * eek
+  }
>  all.equal(oout_ceesigma$hessian[is.bee_or_cee, is.nu_or_sigma], foompter)
[1] TRUE
>  # equation (4.35f) in Aster Theory book
>  foompter <- matrix(0, length(sigma), length(sigma))
>  moobb <- moo$hessian[is.bee.short, is.bee.short]
>  inv <- solve(zwz %*% diag(dee) + diag(length(dee)))
>  for (j in seq_along(sigma)) {
+      eej <- as.numeric(dee_idx == j)
+      for (k in seq_along(sigma)) {
+          eek <- as.numeric(dee_idx == k)
+          foompter[j, k] <-
+              drop(rbind(eej * cee) %*% moobb %*% cbind(eek * cee)) +
+              sum(diag(inv %*% zwz %*% diag(eej) %*% diag(eek))) -
+              2.0 * sum(diag(inv %*% zwz %*% diag(eej * aaa) %*%
+                       inv %*% zwz %*% diag(eek * aaa)))
+      }
+  }
>  all.equal(oout_ceesigma$hessian[is.nu_or_sigma, is.nu_or_sigma], foompter)
[1] TRUE
>  # everything checks
> 
> 
> proc.time()
   user  system elapsed 
   2.04    0.28    2.31