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)
> 
> data(radish)
> 
> pred <- c(0,1,2)
> fam <- c(1,3,2)
> 
> ### try reaster
> 
> rout <- reaster(resp ~ varb + fit : (Site * Region),
+     list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop),
+     pred, fam, varb, id, root, data = radish)
> srout <- summary(rout)
> names(rout)
 [1] "obj"        "fixed"      "random"     "dropped"    "sigma"     
 [6] "nu"         "c"          "b"          "alpha"      "zwz"       
[11] "response"   "origin"     "iterations" "counts"     "deviance"  
[16] "formula"    "call"      
> 
> foo <- new.env(parent = emptyenv())
> bar <- suppressWarnings(try(load("quickle.rda", foo), silent = TRUE))
> if (inherits(bar, "try-error")) {
+     save(srout, file = "quickle.rda")
+ } else {
+     srout$object$iterations <- NULL
+     foo$srout$object$iterations <- NULL
+     print(all.equal(srout, foo$srout, tol = 1e-4))
+ }
[1] TRUE
> 
> alpha.mle <- rout$alpha
> bee.mle <- rout$b
> nu.mle <- rout$sigma^2
> zwz.mle <- rout$zwz
> obj <- rout$obj
> fixed <- rout$fixed
> random <- rout$random
> alphanu.mle <- c(alpha.mle, nu.mle)
> 
> qout <- quickle(alphanu.mle, bee.mle, fixed, random, obj,
+     zwz = zwz.mle, deriv = 2)
> 
> objfun <- function(alphanu) quickle(alphanu, bee.mle, fixed, random,
+     obj, zwz = zwz.mle)$value
> gradfun <- function(alphanu) quickle(alphanu, bee.mle, fixed, random,
+     obj, zwz = zwz.mle, deriv = 1)$gradient
> oout <- optim(alphanu.mle, objfun, gradfun, method = "BFGS", hessian = TRUE)
> all.equal(qout$hessian, oout$hessian, check.attributes = FALSE,
+     tolerance = 0.002)
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   2.00    0.25    2.21