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
bogus.R
library(aster)
data(echinacea)
# cut and paste from help(predict.aster)
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 <- grepl("hdct", as.character(redata$varb))
redata <- data.frame(redata, hdct = as.integer(hdct))
level <- gsub("[0-9]", "", as.character(redata$varb))
redata <- data.frame(redata, level = as.factor(level))
aout <- aster(resp ~ varb + level : (nsloc + ewloc) + hdct : pop,
pred, fam, varb, id, root, data = redata)
# end of cut and paste from help(predict.aster)
nind <- length(unique(redata$id))
nnode <- length(pred)
ncoef <- length(aout$coefficients)
nind
nnode
ncoef
mu <- predict(aout)
mu <- matrix(mu, nrow = nind, ncol = nnode)
my.xi <- matrix(NaN, nrow = nind, ncol = nnode)
for (j in 1:nnode)
if (pred[j] == 0) {
my.xi[ , j] <- mu[ , j]
} else {
my.xi[ , j] <- mu[ , j] / mu[ , pred[j]]
}
old.xi <- predict(aout, model.type = "conditional")
old.xi <- matrix(old.xi, nrow = nind, ncol = nnode)
max(abs(my.xi - old.xi)) > 1.0
# force use of method predict.aster rather than predict.aster.formula
x.bogo <- matrix(1.0, nrow = nind, ncol = nnode)
tricky.xi <- predict.aster(aout, x.bogo, x.bogo, aout$modmat,
model.type = "conditional")
tricky.xi <- matrix(tricky.xi, nrow = nind, ncol = nnode)
all.equal(tricky.xi, my.xi)
### now use new feature !!!
new.xi <- predict(aout, model.type = "conditional", is.always = TRUE)
new.xi <- matrix(new.xi, nrow = nind, ncol = nnode)
all.equal(new.xi, my.xi)
# check gradient
xpred <- old.xi / new.xi
my.xpred <- matrix(NaN, nrow = nind, ncol = nnode)
for (j in 1:nnode)
if (pred[j] == 0) {
my.xpred[ , j] <- 1
} else {
my.xpred[ , j] <- aout$x[ , pred[j]]
}
all.equal(xpred, my.xpred)
pout.new <- predict(aout, model.type = "conditional", gradient = TRUE,
is.always = TRUE)
pout.old <- predict(aout, model.type = "conditional", gradient = TRUE)
dim(pout.new$gradient)
dim(pout.old$gradient)
length(xpred)
my.gradient.old <- sweep(pout.new$gradient, 1, as.vector(xpred), "*")
all.equal(my.gradient.old, pout.old$gradient)
##### check conditional aster model #####
aout <- aster(resp ~ varb + level : (nsloc + ewloc + pop),
pred, fam, varb, id, root, data = redata, type = "conditional")
old.xi <- predict(aout, model.type = "conditional")
old.xi <- matrix(old.xi, nrow = nind, ncol = nnode)
new.xi <- predict(aout, model.type = "conditional", is.always = TRUE)
new.xi <- matrix(new.xi, nrow = nind, ncol = nnode)
all.equal(xpred * new.xi, old.xi)
pout.new <- predict(aout, model.type = "conditional", gradient = TRUE,
is.always = TRUE)
pout.old <- predict(aout, model.type = "conditional", gradient = TRUE)
my.gradient.old <- sweep(pout.new$gradient, 1, as.vector(xpred), "*")
all.equal(my.gradient.old, pout.old$gradient)