tests/valid.R

In aster2: Aster Models

 library(aster2)

 set.seed(42)

 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)
 names(fred)
 names(fred$redata)

 theta <- rnorm(length(fred))
 try(validtheta(fred, theta))
 try(validtheta(fred, theta, mod = "cond"))

 varb <- as.character(fred$redata[[fred$varb.name]]) 
 pred <- fred$repred
 resp <- fred$redata[[fred$response.name]]
 init <- fred$initial
 ypred <- ifelse(pred == 0, init, c(NA, resp)[pred + 1])

 theta.cond <- theta
 theta.unco <- theta
 theta.unco[varb == "n2"] <- (- abs(theta.unco[varb == "n2"]))
 theta.cond[varb == "n2" & ypred > 0] <-
     (- abs(theta.unco[varb == "n2" & ypred > 0]))

 is.validtheta(fred, theta.unco)
 is.validtheta(fred, theta.cond, mod = "cond")
 
 xi <- rnorm(length(fred))
 try(validxi(fred, xi))
 try(validxi(fred, xi, mod = "cond"))

 fixup <- function(xi, varb, ypred) {
     stopifnot(length(xi) == length(varb))
     stopifnot(is.numeric(xi))
     stopifnot(is.finite(xi))
     stopifnot(is.character(varb))
     if (missing(ypred)) {
         ypred <- rep(1, length(xi))
     } else {
         stopifnot(is.numeric(ypred))
         stopifnot(length(ypred) == length(xi))
     }
     todo <- ypred != 0
     xi[todo & varb == "z1"] <- 1 + abs(xi[todo & varb == "z1"])
     xi[todo & varb == "p1"] <- abs(xi[todo & varb == "p1"])
     xi[todo & varb == "b1"] <- 1 / (1 + abs(xi[todo & varb == "b1"]))
     xi[todo & varb == "n2"] <- 1 + xi[todo & varb == "n1"]^2
     my.xi.m1 <- exp(xi[todo & varb == "m1"])
     my.xi.m2 <- exp(xi[todo & varb == "m2"])
     my.xi.m3 <- exp(xi[todo & varb == "m3"])
     xi[todo & varb == "m1"] <- my.xi.m1 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
     xi[todo & varb == "m2"] <- my.xi.m2 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
     xi[todo & varb == "m3"] <- my.xi.m3 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
     return(xi)
 }

 xi.unco <- fixup(xi, varb)
 xi.cond <- fixup(xi, varb, ypred)

 is.validxi(fred, xi.unco)
 is.validxi(fred, xi.cond, mod = "cond")

 ########## now with direction of recession ##########

 delta <- fred$redelta
 all(delta == 0)
 delta[grepl("[mbp]", varb) & resp == 0 & runif(length(fred)) < 1 / 5] <- -1
 fred$redelta <- delta
 is.validasterdata(fred)
 sum(delta != 0)
 
 try(validtheta(fred, theta))
 try(validtheta(fred, theta, mod = "cond"))

 is.validtheta(fred, theta.unco)
 is.validtheta(fred, theta.cond, mod = "cond")
 
 my.sally <- 2 * as.numeric(delta < 0)
 # when loop done
 # my.sally == 2 means delta < 0 and predecessor NOT a. s. zero
 # my.sally == 1 means predecessor a. s. zero (regardless of delta)
 repeat {
     save.my.sally <- my.sally
     foom <- c(0, my.sally)[pred + 1]
     my.sally[foom != 0] <- 1
     if (identical(my.sally, save.my.sally)) break
 }
 my.sally <- my.sally == 1

 try(validxi(fred, xi))
 try(validxi(fred, xi, mod = "cond"))

 fixup <- function(xi, varb, ypred, delta) {
     stopifnot(is.character(varb))
     stopifnot(length(xi) == length(varb))
     stopifnot(is.numeric(xi))
     stopifnot(is.finite(xi))
     stopifnot(is.numeric(ypred))
     stopifnot(is.finite(ypred))
     stopifnot(length(ypred) == length(varb))
     stopifnot(is.numeric(delta))
     stopifnot(is.finite(delta))
     stopifnot(length(delta) == length(varb))
     todo <- ypred != 0
     lower <- delta < 0
     xi[todo & varb == "z1"] <- 1 + abs(xi[todo & varb == "z1"])
     xi[todo & varb == "z1" & lower] <- 1
     xi[todo & varb == "p1"] <- abs(xi[todo & varb == "p1"])
     xi[todo & varb == "p1" & lower] <- 0
     xi[todo & varb == "b1"] <- 1 / (1 + abs(xi[todo & varb == "b1"]))
     xi[todo & varb == "b1" & lower] <- 0
     xi[todo & varb == "n2"] <- 1 + xi[todo & varb == "n1"]^2
     my.xi <- ifelse(lower, 0, exp(xi))
     my.xi.m1 <- my.xi[todo & varb == "m1"]
     my.xi.m2 <- my.xi[todo & varb == "m2"]
     my.xi.m3 <- my.xi[todo & varb == "m3"]
     xi[todo & varb == "m1"] <- my.xi.m1 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
     xi[todo & varb == "m2"] <- my.xi.m2 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
     xi[todo & varb == "m3"] <- my.xi.m3 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
     return(xi)
 }

 xi.unco <- fixup(xi, varb, as.numeric(! my.sally), delta)
 xi.cond <- fixup(xi, varb, ypred, delta)

 is.validxi(fred, xi.unco)
 is.validxi(fred, xi.cond, mod = "cond")