tests/constancy.R

In aster2: Aster Models

 library(aster2)

 ##### parm.type = "theta"

 data(echinacea)
 cmat <- constancy(echinacea, parm.type = "theta")
 nrow(cmat) == 0

 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)

 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))
 
 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)

 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))

 ##### is.same parm.type = "theta"

 theta1 <- rnorm(length(fred$repred))
 theta2 <- rnorm(length(fred$repred))
 ! is.same(theta1, theta2, fred)
 theta2 <- theta1 + t(cmat) %*% rnorm(nrow(cmat))
 theta2 <- as.numeric(theta2)
 is.same(theta1, theta2, fred)

 ##### parm.type = "phi"

 rm(list = ls())

 data(echinacea)
 cmat <- constancy(echinacea, parm.type = "phi")
 nrow(cmat) == 0

 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)

 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))

 # (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)

 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))

 ##### is.same parm.type = "phi"

 phi1 <- rnorm(length(fred$repred))
 phi2 <- rnorm(length(fred$repred))
 ! is.same(phi1, phi2, fred, parm.type = "phi")
 phi2 <- phi1 + t(cmat) %*% rnorm(nrow(cmat))
 phi2 <- as.numeric(phi2)
 is.same(phi1, phi2, fred, parm.type = "phi")