Running 'asterdata.R' [1s]
  Comparing 'asterdata.Rout' to 'asterdata.Rout.save' ... OK
  Running 'constancy.R' [169s]
  Running 'echinacea-cond.R' [2s]
  Comparing 'echinacea-cond.Rout' to 'echinacea-cond.Rout.save' ...66,69d65
< Registered S3 methods overwritten by 'aster2':
<   method              from 
<   as.character.astfam aster
<   print.astfam        aster
  Running 'echinacea.R' [167s]
  Running 'fam-degen.R' [2s]
  Comparing 'fam-degen.Rout' to 'fam-degen.Rout.save' ... OK
  Running 'families.R' [160s]
  Running 'logl.R' [1s]
  Comparing 'logl.Rout' to 'logl.Rout.save' ... OK
  Running 'subset.R' [171s]
  Running 'transform.R' [2s]
  Comparing 'transform.Rout' to 'transform.Rout.save' ... OK
  Running 'valid.R' [1s]
  Comparing 'valid.Rout' to 'valid.Rout.save' ... OK
  Running 'zeros.R' [1s]
  Comparing 'zeros.Rout' to 'zeros.Rout.save' ... OK
Running the tests in 'tests/constancy.R' failed.
Complete output:
  > 
  >  library(aster2)
  Loading required package: Matrix
  This is beta software.
  Unless you need to do aster models with dependence groups,
      use package "aster" instead.
  See help(aster2-package) for differences from package "aster"
      and examples.
  
  > 
  >  ##### parm.type = "theta"
  > 
  >  data(echinacea)
  >  cmat <- constancy(echinacea, parm.type = "theta")
  >  nrow(cmat) == 0
  [1] TRUE
  > 
  >  data(test1)
  >  fred <- asterdata(test1,
  +      vars = c("m1", "m2", "m3", "n1", "n2", "b1", "p1", "z1"),
  +      pred = c(0, 0, 0, 1, 1, 2, 3, 6), group = c(0, 1, 2, 0, 4, 0, 0, 0),
  +      code = c(1, 1, 1, 2, 2, 3, 4, 5),
  +      families = list(fam.multinomial(3), "normal.location.scale",
  +      "bernoulli", "poisson", "zero.truncated.poisson"))
  >  is.validasterdata(fred)
  [1] TRUE
  > 
  >  cmat <- constancy(fred, parm.type = "theta")
  >  multinomial.ind <- sort(unique(fred$redata$id[fred$recode == 1]))
  >  mycmat <- matrix(NA, nrow = length(multinomial.ind), ncol = nrow(fred$redata))
  >  for (ind in seq(along = multinomial.ind))
  +      mycmat[ind, ] <- as.numeric(fred$redata$id == ind & fred$recode == 1)
  >  mycmat <- mycmat[rev(multinomial.ind), ]
  >  identical(cmat, Matrix(mycmat))
  [1] TRUE
  >  
  >  sally <- fred
  >  y <- sally$redata[[sally$response.name]]
  >  v <- as.character(sally$redata$varb)
  >  ind <- sally$redata$id
  >  ind.save <- numeric(0)
  >  cmat.extra <- NULL
  >  is.zero <- rep(FALSE, length(y))
  > 
  >  i.buddy <- sally$regroup
  >  i.buddy[i.buddy == 0] <- NA
  >  i.buddy.buddy <- sally$regroup[i.buddy]
  >  i.buddy.buddy[i.buddy.buddy == 0] <- NA
  >  y.buddy <- y[i.buddy]
  >  y.buddy.buddy <- y[i.buddy.buddy]
  > 
  >  # find an individual for whom to constrain m1 == 0
  >  ok <- y == 0 & v == "m1"
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m2 == 0
  >  ok <- y == 0 & v == "m2" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m3 == 0
  >  ok <- y == 0 & v == "m3" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m1 == 0 and m2 == 0
  >  ok <- y != 0 & y.buddy == 0 & y.buddy.buddy == 0 & (! (ind %in% ind.save))
  >  ok[is.na(ok)] <- FALSE
  >  if (any(ok)) {
  +      k <- seq(along = ok)[ok]
  +      k <- k[1]
  +      j <- i.buddy[k]
  +      i <- i.buddy.buddy[k]
  +      sally$redelta[i] <- (-1)
  +      sally$redelta[j] <- (-1)
  +      is.zero[i] <- TRUE
  +      is.zero[j] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +      foo <- rep(0, length(y))
  +      foo[j] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m1 == 0 and m3 == 0
  >  ok <- y == 0 & y.buddy != 0 & y.buddy.buddy == 0 & (! (ind %in% ind.save))
  >  ok[is.na(ok)] <- FALSE
  >  if (any(ok)) {
  +      k <- seq(along = ok)[ok]
  +      k <- k[1]
  +      j <- i.buddy[k]
  +      i <- i.buddy.buddy[k]
  +      sally$redelta[i] <- (-1)
  +      sally$redelta[k] <- (-1)
  +      is.zero[i] <- TRUE
  +      is.zero[k] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +      foo <- rep(0, length(y))
  +      foo[k] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m2 == 0 and m3 == 0
  >  ok <- y == 0 & y.buddy == 0 & y.buddy.buddy != 0 & (! (ind %in% ind.save))
  >  ok[is.na(ok)] <- FALSE
  >  if (any(ok)) {
  +      k <- seq(along = ok)[ok]
  +      k <- k[1]
  +      j <- i.buddy[k]
  +      i <- i.buddy.buddy[k]
  +      sally$redelta[j] <- (-1)
  +      sally$redelta[k] <- (-1)
  +      is.zero[j] <- TRUE
  +      is.zero[k] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[j] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +      foo <- rep(0, length(y))
  +      foo[k] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain b1 == 0
  >  ok <- y == 0 & v == "b1" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain b1 == 1
  >  ok <- y == 1 & v == "b1" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (+1)
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain p1 == 0
  >  ok <- y == 0 & v == "p1" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain z1 == 1
  >  ok <- y == 1 & v == "z1" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  is.validasterdata(sally)
  [1] TRUE
  > 
  >  i.pred <- sally$repred
  >  i.pred[i.pred == 0] <- NA
  >  anc.is.zero <- rep(FALSE, length(y))
  >  repeat {
  +      anc.is.zero.save <- anc.is.zero
  +      anc.is.zero <- is.zero[i.pred] | anc.is.zero[i.pred]
  +      anc.is.zero[is.na(anc.is.zero)] <- FALSE
  +      if (identical(anc.is.zero, anc.is.zero.save)) break
  +  }
  >  cmat.extra.extra <- matrix(0, sum(anc.is.zero), length(y))
  >  j <- seq(along = anc.is.zero)[anc.is.zero]
  >  for (i in 1:nrow(cmat.extra.extra)) {
  +      cmat.extra.extra[i, j[i]] <- 1
  +  }
  >  
  >  dimnames(cmat.extra) <- NULL
  >  cmat.sally <- constancy(sally, parm.type = "theta")
  >  mycmat.sally <- rbind(cmat, cmat.extra, cmat.extra.extra)
  >  cmat.sally.char <- apply(cmat.sally, 1, paste, collapse = "*")
  >  mycmat.sally.char <- apply(mycmat.sally, 1, paste, collapse = "*")
  >  identical(sort(cmat.sally.char), sort(mycmat.sally.char))
  [1] TRUE
  > 
  >  ##### is.same parm.type = "theta"
  > 
  >  theta1 <- rnorm(length(fred$repred))
  >  theta2 <- rnorm(length(fred$repred))
  >  ! is.same(theta1, theta2, fred)
  [1] TRUE
  >  theta2 <- theta1 + t(cmat) %*% rnorm(nrow(cmat))
  >  theta2 <- as.numeric(theta2)
  >  is.same(theta1, theta2, fred)
  [1] TRUE
  > 
  >  ##### parm.type = "phi"
  > 
  >  rm(list = ls())
  > 
  >  data(echinacea)
  >  cmat <- constancy(echinacea, parm.type = "phi")
  >  nrow(cmat) == 0
  [1] TRUE
  > 
  >  data(test1)
  >  fred <- asterdata(test1,
  +      vars = c("m1", "m2", "m3", "n1", "n2", "b1", "p1", "z1"),
  +      pred = c(0, 0, 0, 1, 1, 2, 3, 6), group = c(0, 1, 2, 0, 4, 0, 0, 0),
  +      code = c(1, 1, 1, 2, 2, 3, 4, 5),
  +      families = list(fam.multinomial(3), "normal.location.scale",
  +      "bernoulli", "poisson", "zero.truncated.poisson"))
  >  is.validasterdata(fred)
  [1] TRUE
  > 
  >  cmat <- constancy(fred, parm.type = "phi")
  >  multinomial.ind <- sort(unique(fred$redata$id[fred$recode == 1]))
  >  mycmat <- matrix(NA, nrow = length(multinomial.ind), ncol = nrow(fred$redata))
  >  for (ind in seq(along = multinomial.ind))
  +      mycmat[ind, ] <- as.numeric(fred$redata$id == ind & fred$recode == 1)
  >  mycmat <- mycmat[rev(multinomial.ind), ]
  >  identical(cmat, Matrix(mycmat))
  [1] TRUE
  > 
  >  # (from the help page for the R function asterdata) "note that the
  >  # conditional canonical parameter vector is always used here [for delta and
  >  # redelta], regardless of whether conditional or unconditional canonical
  >  # affine submodels are to be used
  > 
  >  sally <- fred
  >  y <- sally$redata[[sally$response.name]]
  >  v <- as.character(sally$redata$varb)
  >  ind <- sally$redata$id
  >  ind.save <- numeric(0)
  >  cmat.extra <- NULL
  >  is.zero <- rep(FALSE, length(y))
  > 
  >  i.buddy <- sally$regroup
  >  i.buddy[i.buddy == 0] <- NA
  >  i.buddy.buddy <- sally$regroup[i.buddy]
  >  i.buddy.buddy[i.buddy.buddy == 0] <- NA
  >  y.buddy <- y[i.buddy]
  >  y.buddy.buddy <- y[i.buddy.buddy]
  > 
  >  i.pred <- sally$repred
  >  i.pred[i.pred == 0] <- NA
  >  y.pred <- y[i.pred]
  >  y.pred[is.na(y.pred)] <- 1
  > 
  >  # find an individual for whom to constrain m1 == 0
  >  ok <- y == 0 & v == "m1"
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m2 == 0
  >  ok <- y == 0 & v == "m2" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m3 == 0
  >  ok <- y == 0 & v == "m3" & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m1 == 0 and m2 == 0
  >  ok <- y != 0 & y.buddy == 0 & y.buddy.buddy == 0 & (! (ind %in% ind.save))
  >  ok[is.na(ok)] <- FALSE
  >  if (any(ok)) {
  +      k <- seq(along = ok)[ok]
  +      k <- k[1]
  +      j <- i.buddy[k]
  +      i <- i.buddy.buddy[k]
  +      sally$redelta[i] <- (-1)
  +      sally$redelta[j] <- (-1)
  +      is.zero[i] <- TRUE
  +      is.zero[j] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +      foo <- rep(0, length(y))
  +      foo[j] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m1 == 0 and m3 == 0
  >  ok <- y == 0 & y.buddy != 0 & y.buddy.buddy == 0 & (! (ind %in% ind.save))
  >  ok[is.na(ok)] <- FALSE
  >  if (any(ok)) {
  +      k <- seq(along = ok)[ok]
  +      k <- k[1]
  +      j <- i.buddy[k]
  +      i <- i.buddy.buddy[k]
  +      sally$redelta[i] <- (-1)
  +      sally$redelta[k] <- (-1)
  +      is.zero[i] <- TRUE
  +      is.zero[k] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +      foo <- rep(0, length(y))
  +      foo[k] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain m2 == 0 and m3 == 0
  >  ok <- y == 0 & y.buddy == 0 & y.buddy.buddy != 0 & (! (ind %in% ind.save))
  >  ok[is.na(ok)] <- FALSE
  >  if (any(ok)) {
  +      k <- seq(along = ok)[ok]
  +      k <- k[1]
  +      j <- i.buddy[k]
  +      i <- i.buddy.buddy[k]
  +      sally$redelta[j] <- (-1)
  +      sally$redelta[k] <- (-1)
  +      is.zero[j] <- TRUE
  +      is.zero[k] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[j] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +      foo <- rep(0, length(y))
  +      foo[k] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain b1 == 0
  >  ok <- y == 0 & v == "b1" & y.pred > 0 & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain b1 == 1
  >  ok <- y == 1 & v == "b1" & y.pred > 0 & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- 1
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      foo[i.pred[i]] <- (-1)
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain p1 == 0
  >  ok <- y == 0 & v == "p1" & y.pred > 0 & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      is.zero[i] <- TRUE
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  # find another individual for whom to constrain z1 == 1
  >  ok <- y == 1 & v == "z1" & y.pred > 0 & (! (ind %in% ind.save))
  >  if (any(ok)) {
  +      i <- seq(along = ok)[ok]
  +      i <- i[1]
  +      sally$redelta[i] <- (-1)
  +      ind.save <- c(ind.save, ind[i])
  +      foo <- rep(0, length(y))
  +      foo[i] <- 1
  +      foo[i.pred[i]] <- (-1)
  +      cmat.extra <- rbind(cmat.extra, foo)
  +  }
  >  is.validasterdata(sally)
  [1] TRUE
  > 
  >  anc.is.zero <- rep(FALSE, length(y))
  >  repeat {
  +      anc.is.zero.save <- anc.is.zero
  +      anc.is.zero <- is.zero[i.pred] | anc.is.zero[i.pred]
  +      anc.is.zero[is.na(anc.is.zero)] <- FALSE
  +      if (identical(anc.is.zero, anc.is.zero.save)) break
  +  }
  >  cmat.extra.extra <- matrix(0, sum(anc.is.zero), length(y))
  >  j <- seq(along = anc.is.zero)[anc.is.zero]
  >  for (i in 1:nrow(cmat.extra.extra)) {
  +      cmat.extra.extra[i, j[i]] <- 1
  +  }
  >  
  >  dimnames(cmat.extra) <- NULL
  >  cmat.sally <- constancy(sally, parm.type = "phi")
  >  mycmat.sally <- rbind(cmat, cmat.extra, cmat.extra.extra)
  >  cmat.sally.char <- apply(cmat.sally, 1, paste, collapse = "*")
  >  mycmat.sally.char <- apply(mycmat.sally, 1, paste, collapse = "*")
  >  identical(sort(cmat.sally.char), sort(mycmat.sally.char))
  [1] TRUE
  > 
  >  ##### is.same parm.type = "phi"
  > 
  >  phi1 <- rnorm(length(fred$repred))
  >  phi2 <- rnorm(length(fred$repred))
  >  ! is.same(phi1, phi2, fred, parm.type = "phi")
  [1] TRUE
  >  phi2 <- phi1 + t(cmat) %*% rnorm(nrow(cmat))
  >  phi2 <- as.numeric(phi2)
  >  is.same(phi1, phi2, fred, parm.type = "phi")
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.03    0.18    2.20 
Running the tests in 'tests/echinacea.R' failed.
Complete output:
  > 
  >  #### copied from help page for function library
  > 
  >  pkg <- "aster"
  >  library(pkg, character.only = TRUE)
  > 
  >  set.seed(42)
  > 
  >  #### copied from help page for function aster in package aster
  > 
  >  data(echinacea)
  >  vars <- c("ld02", "ld03", "ld04", "fl02", "fl03", "fl04",
  +      "hdct02", "hdct03", "hdct04")
  >  redata <- reshape(echinacea, varying = list(vars), direction = "long",
  +      timevar = "varb", times = as.factor(vars), v.names = "resp")
  >  redata <- data.frame(redata, root = 1)
  >  pred <- c(0, 1, 2, 1, 2, 3, 4, 5, 6)
  >  fam <- c(1, 1, 1, 1, 1, 1, 3, 3, 3)
  >  hdct <- grep("hdct", as.character(redata$varb))
  >  hdct <- is.element(seq(along = redata$varb), hdct)
  >  redata <- data.frame(redata, hdct = as.integer(hdct))
  >  aout4 <- aster(resp ~ varb + nsloc + ewloc + pop * hdct - pop,
  +      pred, fam, varb, id, root, data = redata)
  >  # summary(aout4)
  > 
  >  beta <- aout4$coefficients
  >  dbeta <- rnorm(length(beta))
  > 
  >  pout <- predict(aout4, parm.type = "canonical", model.type = "conditional",
  +      se.fit = TRUE)
  >  theta <- pout$fit
  >  dtheta <- as.vector(pout$gradient %*% dbeta)
  >  jack.theta <- pout$gradient
  > 
  >  pout <- predict(aout4, parm.type = "canonical", model.type = "unconditional",
  +      se.fit = TRUE)
  >  phi <- pout$fit
  >  dphi <- as.vector(pout$gradient %*% dbeta)
  >  jack.phi <- pout$gradient
  > 
  >  pout <- predict(aout4, parm.type = "mean.value", model.type = "unconditional",
  +      se.fit = TRUE)
  >  mu <- pout$fit
  >  dmu <- as.vector(pout$gradient %*% dbeta)
  >  jack.mu <- pout$gradient
  > 
  >  phony <- matrix(1, nrow = nrow(aout4$x), ncol = ncol(aout4$x))
  >  pout <- predict.aster(aout4, x = phony, root = aout4$root,
  +      modmat = aout4$modmat, parm.type = "mean.value",
  +      model.type = "conditional", se.fit = TRUE)
  >  xi <- pout$fit
  >  dxi <- as.vector(pout$gradient %*% dbeta)
  >  jack.xi <- pout$gradient
  > 
  >  offset <- as.vector(aout4$origin)
  >  modmat <- matrix(aout4$modmat, ncol = length(beta))
  >  tau <- as.numeric(t(modmat) %*% mu)
  >  dtau <- as.numeric(t(modmat) %*% dmu)
  > 
  >  modmat.orig <- model.matrix(resp ~ varb + nsloc + ewloc + pop * hdct - pop,
  +      redata)
  >  beta.orig <- rep(0, ncol(modmat.orig))
  >  beta.orig[match(names(beta), colnames(modmat.orig))] <- beta
  >  names(beta.orig) <- colnames(modmat.orig)
  >  tau.orig <- as.numeric(t(modmat.orig) %*% mu)
  > 
  >  #### copied from help page for function library
  > 
  >  detach(pos = match(paste("package", pkg, sep=":"), search()))
  > 
  >  #### end of stuff from old aster package
  > 
  >  rm(list = setdiff(ls(), c("beta", "theta", "phi", "xi", "mu", "tau",
  +      "offset", "modmat", "dbeta", "dtheta", "dphi", "dxi", "dmu", "dtau",
  +      "jack.theta", "jack.phi", "jack.xi", "jack.mu",
  +      "modmat.orig", "beta.orig", "tau.orig")))
  > 
  >  library(aster2)
  Loading required package: Matrix
  Registered S3 methods overwritten by 'aster2':
    method              from 
    as.character.astfam aster
    print.astfam        aster
  This is beta software.
  Unless you need to do aster models with dependence groups,
      use package "aster" instead.
  See help(aster2-package) for differences from package "aster"
      and examples.
  
  > 
  >  data(echinacea)
  > 
  >  #### saturated
  > 
  >  myphi <- transformSaturated(theta, echinacea, from = "theta", to = "phi")
  >  all.equal(phi, myphi)
  [1] TRUE
  > 
  >  mytheta <- transformSaturated(phi, echinacea, from = "phi", to = "theta")
  >  all.equal(theta, mytheta)
  [1] TRUE
  > 
  >  myxi <- transformSaturated(theta, echinacea, from = "theta", to = "xi")
  >  all.equal(xi, myxi)
  [1] TRUE
  > 
  >  mymu <- transformSaturated(xi, echinacea, from = "xi", to = "mu")
  >  all.equal(mu, mymu)
  [1] TRUE
  > 
  >  #### unconditional from == "beta"
  > 
  >  phi.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "phi", offset = offset)
  >  all.equal(phi, phi.foo)
  [1] TRUE
  > 
  >  theta.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "theta", offset = offset)
  >  all.equal(theta, theta.foo)
  [1] TRUE
  > 
  >  xi.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "xi", offset = offset)
  >  all.equal(xi, xi.foo)
  [1] TRUE
  > 
  >  mu.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "mu", offset = offset)
  >  all.equal(mu, mu.foo)
  [1] TRUE
  > 
  >  tau.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "tau", offset = offset)
  >  all.equal(tau, tau.foo)
  [1] TRUE
  > 
  >  #### unconditional from == "beta" (differential)
  > 
  >  dphi.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "phi", offset = offset, differential = dbeta)
  >  all.equal(dphi, dphi.foo)
  [1] TRUE
  > 
  >  dtheta.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "theta", offset = offset, differential = dbeta)
  >  all.equal(dtheta, dtheta.foo)
  [1] TRUE
  > 
  >  dxi.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "xi", offset = offset, differential = dbeta)
  >  all.equal(dxi, dxi.foo)
  [1] TRUE
  > 
  >  dmu.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "mu", offset = offset, differential = dbeta)
  >  all.equal(dmu, dmu.foo)
  [1] TRUE
  > 
  >  dtau.foo <- transformUnconditional(beta, modmat, echinacea,
  +      from = "beta", to = "tau", offset = offset, differential = dbeta)
  >  all.equal(dtau, dtau.foo)
  [1] TRUE
  > 
  >  #### Jacobian matrices
  > 
  >  my.jack.theta <- jacobian(beta, echinacea, from = "beta", to = "theta",
  +      offset = offset, modmat = modmat, transform = "unconditional")
  >  all.equal(jack.theta, my.jack.theta)
  [1] TRUE
  > 
  >  my.jack.phi <- jacobian(beta, echinacea, from = "beta", to = "phi",
  +      offset = offset, modmat = modmat, transform = "unconditional")
  >  all.equal(jack.phi, my.jack.phi)
  [1] TRUE
  > 
  >  my.jack.xi <- jacobian(beta, echinacea, from = "beta", to = "xi",
  +      offset = offset, modmat = modmat, transform = "unconditional")
  >  all.equal(jack.xi, my.jack.xi)
  [1] TRUE
  > 
  >  my.jack.mu <- jacobian(beta, echinacea, from = "beta", to = "mu",
  +      offset = offset, modmat = modmat, transform = "unconditional")
  >  all.equal(jack.mu, my.jack.mu)
  [1] TRUE
  > 
  >  my.jack.tau <- jacobian(beta, echinacea, from = "beta", to = "tau",
  +      offset = offset, modmat = modmat, transform = "unconditional")
  >  all.equal(t(modmat) %*% jack.mu, my.jack.tau)
  [1] TRUE
  > 
  >  #### unconditional from == "tau"
  > 
  >  beta.foo <- transformUnconditional(tau, modmat, echinacea,
  +      from = "tau", to = "beta", offset = offset)
  >  all.equal(as.vector(beta), beta.foo)
  [1] TRUE
  > 
  >  xi.foo <- transformUnconditional(tau, modmat, echinacea,
  +      from = "tau", to = "xi", offset = offset)
  >  all.equal(as.vector(xi), xi.foo)
  [1] TRUE
  > 
  >  beta.orig.foo <- transformUnconditional(tau.orig, modmat.orig, echinacea,
  +      from = "tau", to = "beta", offset = offset)
  >  delta.orig <- beta.orig.foo - beta.orig
  >  eta.orig <- as.numeric(modmat.orig %*% delta.orig)
  >  consmat <- constancy(echinacea, parm.type = "phi")
  >  consmat.qr <- qr(t(consmat))
  >  eta.orig.resid <- qr.resid(consmat.qr, eta.orig)
  >  all.equal(max(abs(eta.orig.resid)), 0)
  [1] TRUE
  > 
  >  #### unconditional from == "tau" (differential)
  > 
  >  dbeta.foo <- transformUnconditional(tau, modmat, echinacea,
  +      from = "tau", to = "beta", offset = offset, differential = dtau)
  >  all.equal(dbeta, dbeta.foo)
  [1] TRUE
  > 
  >  dxi.foo <- transformUnconditional(tau, modmat, echinacea,
  +      from = "tau", to = "xi", offset = offset, differential = dtau)
  >  all.equal(as.vector(dxi), dxi.foo)
  [1] TRUE
  > 
  >  dbeta.orig <- rnorm(length(tau.orig))
  >  dtau.orig <- transformUnconditional(beta.orig, modmat.orig, echinacea,
  +      from = "beta", to = "tau", offset = offset, differential = dbeta.orig)
  >  dbeta.orig.foo <- transformUnconditional(tau.orig, modmat.orig, echinacea,
  +      from = "tau", to = "beta", offset = offset, differential = dtau.orig)
  >  # since betas are meaningless, cannot compare them
  >  dphi.orig <- as.numeric(modmat.orig %*% dbeta.orig)
  >  dphi.orig.foo <- as.numeric(modmat.orig %*% dbeta.orig.foo)
  >  is.same(dphi.orig, dphi.orig.foo, echinacea, parm.type = "phi")
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
    11.34    1.71   13.18 
Running the tests in 'tests/families.R' failed.
Complete output:
  > 
  >  library(aster2)
  Loading required package: Matrix
  This is beta software.
  Unless you need to do aster models with dependence groups,
      use package "aster" instead.
  See help(aster2-package) for differences from package "aster"
      and examples.
  
  > 
  >  set.seed(42)
  > 
  >  epsilon <- 1e-6
  > 
  >  ##### Bernoulli #####
  > 
  >  theta <- seq(-3, 3, 0.1)
  > 
  >  cumfun <- function(theta)
  +      ifelse(theta > 0, theta + log1p(exp(- theta)), log1p(exp(theta)))
  > 
  >  moofun <- function(theta)
  +      ifelse(theta > 0, 1 / (1 + exp(- theta)), exp(theta) / (1 + exp(theta)))
  > 
  >  voofun <- function(theta) {
  +      pp <- moofun(theta)
  +      qq <- moofun(- theta)
  +      pp * qq
  +  }
  > 
  >  thoofun <- function(theta) {
  +      pp <- moofun(theta)
  +      qq <- moofun(- theta)
  +      - pp * qq * tanh(theta / 2)
  +  }
  > 
  >  zeroth <- mapply(function(theta) cumulant(theta, fam.bernoulli())$zeroth,
  +      theta)
  > 
  >  all.equal(zeroth, cumfun(theta))
  [1] TRUE
  > 
  >  first <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
  +      deriv = 1)$first, theta)
  > 
  >  all.equal(first, moofun(theta))
  [1] TRUE
  > 
  >  second <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
  +      deriv = 2)$second, theta)
  > 
  >  all.equal(second, voofun(theta))
  [1] TRUE
  > 
  >  third <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
  +      deriv = 3)$third, theta)
  > 
  >  all.equal(third, thoofun(theta))
  [1] TRUE
  > 
  >  grad <- moofun(theta)
  >  fdgrad <- (cumfun(theta + epsilon) - cumfun(theta - epsilon)) / (2 * epsilon)
  >  all.equal(grad, fdgrad)
  [1] TRUE
  > 
  >  hess <- voofun(theta)
  >  fdhess <- (moofun(theta + epsilon) - moofun(theta - epsilon)) / (2 * epsilon)
  >  all.equal(hess, fdhess)
  [1] TRUE
  > 
  >  thss <- thoofun(theta)
  >  fdth <- (voofun(theta + epsilon) - voofun(theta - epsilon)) / (2 * epsilon)
  >  all.equal(thss, fdth)
  [1] TRUE
  > 
  >  foo <- mapply(function(xi) link(xi, fam.bernoulli(), deriv = 0)$zeroth, first)
  > 
  >  all.equal(theta, foo)
  [1] TRUE
  > 
  >  foo <- mapply(function(xi) link(xi, fam.bernoulli(), deriv = 1)$first, first)
  > 
  >  all.equal(second, 1 / foo)
  [1] TRUE
  > 
  >  ##### Poisson #####
  > 
  >  zeroth <- mapply(function(theta) cumulant(theta, fam.poisson())$zeroth, theta)
  > 
  >  all.equal(zeroth, exp(theta))
  [1] TRUE
  > 
  >  first <- mapply(function(theta) cumulant(theta, fam.poisson(),
  +      deriv = 1)$first, theta)
  > 
  >  all.equal(first, exp(theta))
  [1] TRUE
  > 
  >  second <- mapply(function(theta) cumulant(theta, fam.poisson(),
  +      deriv = 2)$second, theta)
  > 
  >  all.equal(second, exp(theta))
  [1] TRUE
  > 
  >  third <- mapply(function(theta) cumulant(theta, fam.poisson(),
  +      deriv = 3)$third, theta)
  > 
  >  all.equal(third, exp(theta))
  [1] TRUE
  > 
  >  foo <- mapply(function(xi) link(xi, fam.poisson(), deriv = 0)$zeroth, first)
  > 
  >  all.equal(theta, foo)
  [1] TRUE
  > 
  >  foo <- mapply(function(xi) link(xi, fam.poisson(), deriv = 1)$first, first)
  > 
  >  all.equal(second, 1 / foo)
  [1] TRUE
  > 
  >  ##### zero-truncated Poisson #####
  > 
  >  cumfun <- function(theta) {
  +      m <- exp(theta)
  +      result <- m + log1p(- exp(- m))
  +      bar <- m / 2 * (1 + m / 3 * (1 + m / 4 * (1 + m / 5 * (1 + m / 6 *
  +          (1 + m / 7 * (1 + m / 8))))))
  +      inies <- theta < (- 4)
  +      result[inies] <- (theta + log1p(bar))[inies]
  +      return(result)
  +  }
  > 
  >  moofun <- function(theta) {
  +      m <- exp(theta)
  +      result <- m / (1 - exp(- m))
  +      bar <- m / 2 * (1 + m / 3 * (1 + m / 4 * (1 + m / 5 * (1 + m / 6 *
  +          (1 + m / 7 * (1 + m / 8))))))
  +      inies <- theta < (- 4)
  +      result[inies] <- (m + 1 / (1 + bar))[inies]
  +      return(result)
  +  }
  > 
  >  voofun <- function(theta) {
  +      m <- exp(theta)
  +      mu <- moofun(theta)
  +      return(mu * (1 + m - mu))
  +  }
  > 
  >  thoofun <- function(theta) {
  +      m <- exp(theta)
  +      mu <- moofun(theta)
  +      return(mu * ((1 + m - mu) * (1 + m - 2 * mu) + m))
  +  }
  > 
  >  theta <- seq(-12, 2, 0.1)
  > 
  >  zeroth <- mapply(function(theta)
  +      cumulant(theta, fam.zero.truncated.poisson())$zeroth, theta)
  > 
  >  all.equal(zeroth, cumfun(theta))
  [1] TRUE
  > 
  >  first <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
  +      deriv = 1)$first, theta)
  > 
  >  all.equal(first, moofun(theta))
  [1] TRUE
  > 
  >  second <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
  +      deriv = 2)$second, theta)
  > 
  >  all.equal(second, voofun(theta))
  [1] TRUE
  > 
  >  third <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
  +      deriv = 3)$third, theta)
  > 
  >  all.equal(third, thoofun(theta))
  [1] TRUE
  > 
  >  grad <- moofun(theta)
  >  fdgrad <- (cumfun(theta + epsilon) - cumfun(theta - epsilon)) / (2 * epsilon)
  >  all.equal(grad, fdgrad)
  [1] TRUE
  > 
  >  hess <- voofun(theta)
  >  fdhess <- (moofun(theta + epsilon) - moofun(theta - epsilon)) / (2 * epsilon)
  >  all.equal(hess, fdhess)
  [1] TRUE
  > 
  >  thss <- thoofun(theta)
  >  fdth <- (voofun(theta + epsilon) - voofun(theta - epsilon)) / (2 * epsilon)
  >  all.equal(thss, fdth)
  [1] TRUE
  > 
  >  foo <- mapply(function(xi) link(xi, fam.zero.truncated.poisson(),
  +      deriv = 0)$zeroth, first)
  > 
  >  all.equal(theta, foo)
  [1] TRUE
  > 
  >  foo <- mapply(function(xi) link(xi, fam.zero.truncated.poisson(),
  +      deriv = 1)$first, first)
  > 
  >  all.equal(second, 1 / foo)
  [1] TRUE
  > 
  >  ##### Multinomial #####
  > 
  >  d <- 4
  > 
  >  cumfun <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      theta.max = max(theta)
  +      j <- seq(along = theta)[theta == theta.max]
  +      j <- j[1]
  +      return(theta.max + log1p(sum(exp(theta - theta.max)[-j])))
  +  }
  > 
  >  moofun <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      theta.max = max(theta)
  +      return(exp(theta - theta.max) / sum(exp(theta - theta.max)))
  +  }
  > 
  >  voofun <- function(theta) {
  +      mu <- moofun(theta)
  +      diag(mu) - outer(mu, mu)
  +  }
  > 
  >  thoofun <- function(theta) {
  +      mu <- moofun(theta)
  +      result <- 2 * outer(mu, outer(mu, mu))
  +      for (i in 1:d) {
  +          for (j in 1:d) {
  +              if (i != j) {
  +                  qux <- mu[i] * mu[j]
  +                  result[i, i, j] <- result[i, i, j] - qux
  +                  result[i, j, i] <- result[i, j, i] - qux
  +                  result[j, i, i] <- result[j, i, i] - qux
  +              }
  +          }
  +          qux <- mu[i]
  +          result[i, i, i] <- qux * (1 - qux) * (1 - 2 * qux)
  +      }
  +      return(result)
  +  }
  > 
  >  theta <- matrix(rnorm(d * 25), ncol = d)
  > 
  >  itself <- apply(theta, 1, cumfun)
  > 
  >  grad <- t(apply(theta, 1, moofun))
  >  all.equal(rowSums(grad), rep(1, nrow(grad)))
  [1] TRUE
  > 
  >  fred <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      result <- double(d)
  +      for (i in 1:d) {
  +          theta.plus <- theta
  +          theta.plus[i] <- theta.plus[i] + epsilon
  +          theta.minus <- theta
  +          theta.minus[i] <- theta.minus[i] - epsilon
  +          result[i] <- (cumfun(theta.plus) - cumfun(theta.minus)) / (2 * epsilon)
  +      }
  +      return(result)
  +  }
  >  fdgrad <- t(apply(theta, 1, fred))
  >  all.equal(grad, fdgrad)
  [1] TRUE
  > 
  >  hess <- apply(theta, 1, voofun)
  >  hess <- array(as.vector(hess), dim = c(d, d, ncol(hess)))
  >  dim(hess)
  [1]  4  4 25
  > 
  >  fred <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      result <- matrix(NA, d, d)
  +      for (i in 1:d) {
  +          theta.plus <- theta
  +          theta.plus[i] <- theta.plus[i] + epsilon
  +          theta.minus <- theta
  +          theta.minus[i] <- theta.minus[i] - epsilon
  +          result[i, ] <- (moofun(theta.plus) - moofun(theta.minus)) /
  +              (2 * epsilon)
  +      }
  +      return(result)
  +  }
  > 
  >  fdhess <- apply(theta, 1, fred)
  >  fdhess <- array(as.vector(fdhess), dim = c(d, d, ncol(fdhess)))
  >  all.equal(hess, fdhess)
  [1] TRUE
  > 
  >  thss <- apply(theta, 1, thoofun)
  >  thss <- array(as.vector(thss), dim = c(d, d, d, ncol(thss)))
  >  dim(thss)
  [1]  4  4  4 25
  > 
  >  fred <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      result <- array(NA, dim = c(d, d, d))
  +      for (i in 1:d) {
  +          theta.plus <- theta
  +          theta.plus[i] <- theta.plus[i] + epsilon
  +          theta.minus <- theta
  +          theta.minus[i] <- theta.minus[i] - epsilon
  +          result[i, , ] <- (voofun(theta.plus) - voofun(theta.minus)) /
  +              (2 * epsilon)
  +      }
  +      return(result)
  +  }
  >  fdthss <- apply(theta, 1, fred)
  >  fdthss <- array(as.vector(fdthss), dim = c(d, d, d, ncol(fdthss)))
  >  all.equal(thss, fdthss)
  [1] TRUE
  > 
  >  cumfun <- function(theta) cumulant(theta, fam.multinomial(d))$zeroth
  >  moofun <- function(theta) cumulant(theta, fam.multinomial(d), deriv = 1)$first
  >  voofun <- function(theta) cumulant(theta, fam.multinomial(d), deriv = 2)$second
  >  thoofun <- function(theta) cumulant(theta, fam.multinomial(d), deriv = 3)$third
  > 
  >  all.equal(itself, apply(theta, 1, cumfun))
  [1] TRUE
  >  all.equal(grad, t(apply(theta, 1, moofun)))
  [1] TRUE
  >  chess <- apply(theta, 1, voofun)
  >  chess <- array(as.vector(chess), dim = c(d, d, ncol(chess)))
  >  all.equal(hess, chess)
  [1] TRUE
  >  cthss <- apply(theta, 1, thoofun)
  >  cthss <- array(as.vector(cthss), dim = c(d, d, d, ncol(cthss)))
  >  all.equal(cthss, thss)
  [1] TRUE
  > 
  >  foo <- apply(grad, 1, function(xi) link(xi, fam.multinomial(d),
  +      deriv = 0)$zeroth)
  >  bar <- apply(t(foo), 1, function(theta) cumulant(theta, fam.multinomial(d),
  +      deriv = 1)$first)
  >  all.equal(t(bar), grad)
  [1] TRUE
  > 
  >  baz <- apply(grad, 1, function(xi) link(xi, fam.multinomial(d),
  +      deriv = 1)$first)
  >  baz <- array(as.vector(baz), dim = c(d, d, nrow(theta)))
  >  fdjackfun <- function(xi) {
  +      result <- matrix(NA, d, d)
  +      for (i in 1:d) {
  +          xi.plus <- xi
  +          xi.plus[i] <- xi[i] + epsilon
  +          xi.minus <- xi
  +          xi.minus[i] <- xi[i] - epsilon
  +          result[ , i] <- (link(xi.plus, fam.multinomial(d))$zeroth -
  +              link(xi.minus, fam.multinomial(d))$zeroth) / (2 * epsilon)
  +      }
  +      return(result)
  +  }
  >  qux <- apply(grad, 1, fdjackfun)
  >  qux <- array(as.vector(qux), dim = c(d, d, nrow(theta)))
  >  all.equal(baz, qux)
  [1] TRUE
  > 
  >  ##### Normal location-scale #####
  > 
  >  d <- 2
  > 
  >  cumfun <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      stopifnot(theta[2] < 0)
  +      - theta[1]^2 / (4 * theta[2]) + (1 / 2) * log(- 1 / (2 * theta[2]))
  +  }
  > 
  >  moofun <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      stopifnot(theta[2] < 0)
  +      r1 <- (- theta[1] / (2 * theta[2]))
  +      r2 <- theta[1]^2 / (4 * theta[2]^2) - 1 / (2 * theta[2])
  +      return(c(r1, r2))
  +  }
  > 
  >  voofun <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      stopifnot(theta[2] < 0)
  +      r11 <- (- 1 / (2 * theta[2]))
  +      r12 <- theta[1] / (2 * theta[2]^2)
  +      r22 <- (- theta[1]^2 / (2 * theta[2]^3)) + 1 / (2 * theta[2]^2)
  +      return(matrix(c(r11, r12, r12, r22), 2, 2))
  +  }
  > 
  >  thoofun <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      stopifnot(theta[2] < 0)
  +      r111 <- 0
  +      r112 <- 1 / (2 * theta[2]^2)
  +      r122 <- (- theta[1] / theta[2]^3)
  +      r222 <- 3 * theta[1]^2 / (2 * theta[2]^4) - 1 / theta[2]^3
  +      return(array(c(r111, r112, r112, r122, r112, r122, r122, r222), dim = c(2, 2, 2)))
  +  }
  > 
  >  theta <- cbind(rnorm(25), - rexp(25))
  > 
  >  itself <- apply(theta, 1, cumfun)
  > 
  >  grad <- t(apply(theta, 1, moofun))
  > 
  >  fred <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      result <- double(d)
  +      for (i in 1:d) {
  +          theta.plus <- theta
  +          theta.plus[i] <- theta.plus[i] + epsilon
  +          theta.minus <- theta
  +          theta.minus[i] <- theta.minus[i] - epsilon
  +          result[i] <- (cumfun(theta.plus) - cumfun(theta.minus)) / (2 * epsilon)
  +      }
  +      return(result)
  +  }
  >  fdgrad <- t(apply(theta, 1, fred))
  >  all.equal(grad, fdgrad)
  [1] TRUE
  > 
  >  hess <- apply(theta, 1, voofun)
  >  hess <- array(as.vector(hess), dim = c(d, d, ncol(hess)))
  >  dim(hess)
  [1]  2  2 25
  > 
  >  fred <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      result <- matrix(NA, d, d)
  +      for (i in 1:d) {
  +          theta.plus <- theta
  +          theta.plus[i] <- theta.plus[i] + epsilon
  +          theta.minus <- theta
  +          theta.minus[i] <- theta.minus[i] - epsilon
  +          result[i, ] <- (moofun(theta.plus) - moofun(theta.minus)) /
  +              (2 * epsilon)
  +      }
  +      return(result)
  +  }
  > 
  >  fdhess <- apply(theta, 1, fred)
  >  fdhess <- array(as.vector(fdhess), dim = c(d, d, ncol(fdhess)))
  >  all.equal(hess, fdhess)
  [1] TRUE
  > 
  >  thss <- apply(theta, 1, thoofun)
  >  thss <- array(as.vector(thss), dim = c(d, d, d, ncol(thss)))
  >  dim(thss)
  [1]  2  2  2 25
  > 
  >  fred <- function(theta) {
  +      stopifnot(length(theta) == d)
  +      result <- array(NA, dim = c(d, d, d))
  +      for (i in 1:d) {
  +          theta.plus <- theta
  +          theta.plus[i] <- theta.plus[i] + epsilon
  +          theta.minus <- theta
  +          theta.minus[i] <- theta.minus[i] - epsilon
  +          result[i, , ] <- (voofun(theta.plus) - voofun(theta.minus)) /
  +              (2 * epsilon)
  +      }
  +      return(result)
  +  }
  > 
  >  fdthss <- apply(theta, 1, fred)
  >  fdthss <- array(as.vector(fdthss), dim = c(d, d, d, ncol(fdthss)))
  >  all.equal(thss, fdthss)
  [1] TRUE
  > 
  >  cumfun <- function(theta) cumulant(theta, fam.normal.location.scale())$zeroth
  >  moofun <- function(theta) cumulant(theta, fam.normal.location.scale(),
  +      deriv = 1)$first
  >  voofun <- function(theta) cumulant(theta, fam.normal.location.scale(),
  +      deriv = 2)$second
  >  thoofun <- function(theta) cumulant(theta, fam.normal.location.scale(),
  +      deriv = 3)$third
  > 
  >  all.equal(itself, apply(theta, 1, cumfun))
  [1] TRUE
  >  all.equal(grad, t(apply(theta, 1, moofun)))
  [1] TRUE
  >  chess <- apply(theta, 1, voofun)
  >  chess <- array(as.vector(chess), dim = c(d, d, ncol(chess)))
  >  all.equal(hess, chess)
  [1] TRUE
  >  cthss <- apply(theta, 1, thoofun)
  >  cthss <- array(as.vector(cthss), dim = c(d, d, d, ncol(cthss)))
  >  all.equal(cthss, thss)
  [1] TRUE
  > 
  >  foo <- apply(grad, 1, function(xi) link(xi, fam.normal.location.scale(),
  +      deriv = 0)$zeroth)
  >  foo <- t(foo)
  >  all.equal(theta, foo)
  [1] TRUE
  > 
  >  foo <- apply(grad, 1, function(xi) link(xi, fam.normal.location.scale(),
  +      deriv = 1)$first)
  >  foo <- array(as.vector(foo), dim = c(2, 2, nrow(theta)))
  >  is.ok <- as.logical(rep(0, nrow(theta)))
  >  for (i in 1:nrow(theta))
  +      is.ok[i] <- all.equal(foo[ , , i] %*% hess[ , , i], diag(2))
  >  all(is.ok)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.00    0.25    2.23 
Running the tests in 'tests/subset.R' failed.
Complete output:
  > 
  >  library(aster2)
  Loading required package: Matrix
  This is beta software.
  Unless you need to do aster models with dependence groups,
      use package "aster" instead.
  See help(aster2-package) for differences from package "aster"
      and examples.
  
  > 
  >  data(echinacea)
  > 
  >  foo <- echinacea$redata$pop
  >  bar <- match(levels(foo), as.character(foo))
  >  baz <- echinacea$redata$id %in% echinacea$redata$id[bar]
  > 
  >  out <- subset(echinacea, baz)
  > 
  >  all(sort(unique(out$redata$id)) == sort(bar))
  [1] TRUE
  > 
  >  opred <- echinacea$repred[baz]
  >  ipred <- seq(along = echinacea$repred)[baz]
  >  all(c(0, ipred)[1 + out$repred] == opred)
  [1] TRUE
  > 
  >  opred <- echinacea$regroup[baz]
  >  all(c(0, ipred)[1 + out$regroup] == opred)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.21    0.10    1.32