R Under development (unstable) (2024-09-15 r87152 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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(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 
   1.43    0.14    1.56