tests/mlogl-c-export.R

In aster: Aster Models

# copied from ../man/aster.Rd examples section

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

# redo with new function

 theta <- predict(aout, parm.type = "canonical", model.type = "conditional")
 phi <- predict(aout, parm.type = "canonical", model.type = "unconditional")

 nind <- nrow(echinacea)
 nnode <- length(pred)

 aster:::setfam(fam.default())

 val.theta <- .C(aster:::C_aster_export_exerciser,
     nind = as.integer(nind),
     nnode = as.integer(nnode),
     pred = as.integer(pred),
     fam = as.integer(fam),
     parm = as.double(theta),
     root = as.double(redata$root),
     response = as.double(redata$resp),
     is.unco = FALSE,
     value = double(1))$value

 val.phi <- .C(aster:::C_aster_export_exerciser,
     nind = as.integer(nind),
     nnode = as.integer(nnode),
     pred = as.integer(pred),
     fam = as.integer(fam),
     parm = as.double(phi),
     root = as.double(redata$root),
     response = as.double(redata$resp),
     is.unco = TRUE,
     value = double(1))$value

 aster:::clearfam()

 all.equal(val.theta, aout$deviance / 2)
 all.equal(val.phi, aout$deviance / 2)