https://github.com/cran/aster
Tip revision: 7016e6b97f24943bdab11323884baf030f38260b authored by Charles J. Geyer on 06 July 2018, 07:20:08 UTC
version 1.0-2
version 1.0-2
Tip revision: 7016e6b
quickle.Rout.save
R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: i686-pc-linux-gnu (32-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
>
> 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.888 0.056 2.952