R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.2.0 beta (2005-09-26 r35681)
ISBN 3-900051-07-0

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 a HTML browser interface to help.
Type 'q()' to quit R.

> 
>  library(aster)
Loading required package: trust
> 
>  data(echinacea)
>  names(echinacea)
 [1] "hdct02" "hdct03" "hdct04" "pop"    "ewloc"  "nsloc"  "ld02"   "fl02"  
 [9] "ld03"   "fl03"   "ld04"   "fl04"  
>  sapply(echinacea, class)
   hdct02    hdct03    hdct04       pop     ewloc     nsloc      ld02      fl02 
"integer" "integer" "integer"  "factor" "integer" "integer" "integer" "integer" 
     ld03      fl03      ld04      fl04 
"integer" "integer" "integer" "integer" 
> 
>  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")
>  names(redata)
[1] "pop"   "ewloc" "nsloc" "varb"  "resp"  "id"   
> 
>  redata <- data.frame(redata, root = 1)
>  names(redata)
[1] "pop"   "ewloc" "nsloc" "varb"  "resp"  "id"    "root" 
> 
>  pred <- c(0, 1, 2, 1, 2, 3, 4, 5, 6)
>  fam <- c(1, 1, 1, 1, 1, 1, 3, 3, 3)
> 
>  hdct <- grep("hdct", as.character(redata$varb))
>  hdct <- is.element(seq(along = redata$varb), hdct)
>  redata <- data.frame(redata, hdct = as.integer(hdct))
>  names(redata)
[1] "pop"   "ewloc" "nsloc" "varb"  "resp"  "id"    "root"  "hdct" 
> 
> ###################################################
> ### chunk number 16: fit-4
> ###################################################
>  aout4 <- aster(resp ~ varb + nsloc + ewloc + pop * hdct - pop,
+      pred, fam, varb, id, root, data = redata)
>  summary(aout4, show.graph = TRUE)

Call:
aster.formula(formula = resp ~ varb + nsloc + ewloc + pop * hdct - 
    pop, pred = pred, fam = fam, varvar = varb, idvar = id, root = root, 
    data = redata)


Graphical Model:
 variable predecessor family          
 ld02     root        bernoulli       
 ld03     ld02        bernoulli       
 ld04     ld03        bernoulli       
 fl02     ld02        bernoulli       
 fl03     ld03        bernoulli       
 fl04     ld04        bernoulli       
 hdct02   fl02        non.zero.poisson
 hdct03   fl03        non.zero.poisson
 hdct04   fl04        non.zero.poisson

                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -1.463638   0.178183  -8.214  < 2e-16 ***
varbfl03        -0.328488   0.265142  -1.239 0.215377    
varbfl04        -0.349519   0.240775  -1.452 0.146601    
varbhdct02       1.808447   0.258014   7.009 2.40e-12 ***
varbhdct03       1.825590   0.210824   8.659  < 2e-16 ***
varbhdct04       2.338955   0.197099  11.867  < 2e-16 ***
varbld02        -0.989146   0.311742  -3.173 0.001509 ** 
varbld03         0.776572   0.392815   1.977 0.048049 *  
varbld04         3.918323   0.330729  11.848  < 2e-16 ***
nsloc            0.013600   0.001729   7.867 3.64e-15 ***
ewloc            0.006473   0.001725   3.753 0.000174 ***
popEriley:hdct  -0.173902   0.089753  -1.938 0.052676 .  
popLf:hdct      -0.157706   0.096580  -1.633 0.102490    
popNessman:hdct -0.315681   0.139822  -2.258 0.023962 *  
popNWLF:hdct    -0.109268   0.083325  -1.311 0.189740    
popSPP:hdct      0.023541   0.086552   0.272 0.785630    
popStevens:hdct -0.127403   0.089518  -1.423 0.154677    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Original predictor variables dropped (aliased)
     hdct 
> 
> 
> ###################################################
> ### chunk number 18: conf-level
> ###################################################
>  conf.level <- 0.95
>  crit <- qnorm((1 + conf.level) / 2)
> 
> 
> ###################################################
> ### chunk number 22: predict-newdata
> ###################################################
>  newdata <- data.frame(pop = levels(echinacea$pop))
>  for (v in vars)
+      newdata[[v]] <- 1
>  newdata$root <- 1
>  newdata$ewloc <- 0
>  newdata$nsloc <- 0
> 
> 
> ###################################################
> ### chunk number 23: predict-newdata-reshape
> ###################################################
>  renewdata <- reshape(newdata, varying = list(vars),
+       direction = "long", timevar = "varb", times = as.factor(vars),
+       v.names = "resp")
>  hdct <- grep("hdct", as.character(renewdata$varb))
>  hdct <- is.element(seq(along = renewdata$varb), hdct)
>  renewdata <- data.frame(renewdata, hdct = as.integer(hdct))
>  names(redata)
[1] "pop"   "ewloc" "nsloc" "varb"  "resp"  "id"    "root"  "hdct" 
>  names(renewdata)
[1] "pop"   "root"  "ewloc" "nsloc" "varb"  "resp"  "id"    "hdct" 
> 
> 
> ###################################################
> ### chunk number 24: make-amat
> ###################################################
>  nind <- nrow(newdata)
>  nnode <- length(vars)
>  amat <- array(0, c(nind, nnode, nind))
>  for (i in 1:nind)
+      amat[i , grep("hdct", vars), i] <- 1
> 
> 
> ###################################################
> ### chunk number 30: tau-4-amat
> ###################################################
>  pout4 <- predict(aout4, varvar = varb, idvar = id, root = root,
+      newdata = renewdata, se.fit = TRUE, amat = amat)
> 
> 
> ###################################################
> ### chunk number 36: make-theta
> ###################################################
>  theta.hat <- predict(aout4, model.type = "cond", parm.type = "canon")
>  theta.hat <- matrix(theta.hat, nrow = nrow(aout4$x), ncol = ncol(aout4$x))
>  fit.hat <- pout4$fit
>  beta.hat <- aout4$coefficients
> 
> 
> ###################################################
> ### chunk number 37: make-root-etc
> ###################################################
>  root <- aout4$root
>  modmat <- aout4$modmat
>  modmat.pred <- pout4$modmat
>  x.pred <- matrix(1, nrow = dim(modmat.pred)[1], ncol = dim(modmat.pred)[2])
>  root.pred <- x.pred
> 
> 
> ###################################################
> ### save and restore
> ###################################################
> 
>  quux <- list(fit.hat = fit.hat, theta.hat = theta.hat,
+      pred.orig = pred, fam.orig = fam, root.orig = root,
+      modmat.orig = modmat, beta.hat = beta.hat,
+      x.pred = x.pred, root.pred = root.pred, modmat.pred = modmat.pred,
+      amat.orig = amat, crit.orig = crit)
>  save(quux, file = "tutor.save.bork.bork.bork.RData")
>  rm(list = ls(all.names = TRUE))
>  load(file = "tutor.save.bork.bork.bork.RData")
>  unlink("tutor.save.bork.bork.bork.RData")
>  attach(quux)
>  ls(pos = 2)
 [1] "amat.orig"   "beta.hat"    "crit.orig"   "fam.orig"    "fit.hat"    
 [6] "modmat.orig" "modmat.pred" "pred.orig"   "root.orig"   "root.pred"  
[11] "theta.hat"   "x.pred"     
> 
> ###################################################
> ### chunk number 38: doit
> ###################################################
> 
>  set.seed(42)
>  nboot <- 5
>  cover <- matrix(0, nboot, length(fit.hat))
>  for (iboot in 1:nboot) {
+      xstar <- raster(theta.hat, pred.orig, fam.orig, root.orig)
+      aout4star <- aster(xstar, root.orig, pred.orig, fam.orig,
+          modmat.orig, beta.hat)
+      pout4star <- predict(aout4star, x.pred, root.pred, modmat.pred,
+          amat.orig, se.fit = TRUE)
+      upper <- pout4star$fit + crit.orig * pout4star$se.fit
+      lower <- pout4star$fit - crit.orig * pout4star$se.fit
+      cover[iboot, ] <- as.numeric(lower <= fit.hat & fit.hat <= upper)
+  }
>  pboot <- apply(cover, 2, mean)
>  pboot.se <- sqrt(pboot * (1 - pboot) / nboot)
>  print(cbind(pboot, pboot.se))
     pboot  pboot.se
[1,]   0.8 0.1788854
[2,]   0.8 0.1788854
[3,]   1.0 0.0000000
[4,]   0.8 0.1788854
[5,]   1.0 0.0000000
[6,]   1.0 0.0000000
[7,]   1.0 0.0000000
> 
> ###################################################
> ### chunk number 38: doit (again)
> ###################################################
> 
>  set.seed(42)
>  nboot <- 5
>  cover <- matrix(0, nboot, length(fit.hat))
>  for (iboot in 1:nboot) {
+      xstar <- raster(theta.hat, pred.orig, fam.orig, root.orig)
+      aout4star <- aster(xstar, root.orig, pred.orig, fam.orig,
+          modmat.orig, beta.hat, method = "CG")
+      pout4star <- predict(aout4star, x.pred, root.pred, modmat.pred,
+          amat.orig, se.fit = TRUE)
+      upper <- pout4star$fit + crit.orig * pout4star$se.fit
+      lower <- pout4star$fit - crit.orig * pout4star$se.fit
+      cover[iboot, ] <- as.numeric(lower <= fit.hat & fit.hat <= upper)
+  }
>  pboot <- apply(cover, 2, mean)
>  pboot.se <- sqrt(pboot * (1 - pboot) / nboot)
>  print(cbind(pboot, pboot.se))
     pboot  pboot.se
[1,]   0.8 0.1788854
[2,]   0.8 0.1788854
[3,]   1.0 0.0000000
[4,]   0.8 0.1788854
[5,]   1.0 0.0000000
[6,]   1.0 0.0000000
[7,]   1.0 0.0000000
> 
> ###################################################
> ### chunk number 38: doit (yet again)
> ###################################################
> 
>  set.seed(42)
>  nboot <- 5
>  cover <- matrix(0, nboot, length(fit.hat))
>  for (iboot in 1:nboot) {
+      xstar <- raster(theta.hat, pred.orig, fam.orig, root.orig)
+      aout4star <- aster(xstar, root.orig, pred.orig, fam.orig,
+          modmat.orig, beta.hat, method = "L-B")
+      pout4star <- predict(aout4star, x.pred, root.pred, modmat.pred,
+          amat.orig, se.fit = TRUE)
+      upper <- pout4star$fit + crit.orig * pout4star$se.fit
+      lower <- pout4star$fit - crit.orig * pout4star$se.fit
+      cover[iboot, ] <- as.numeric(lower <= fit.hat & fit.hat <= upper)
+  }
>  pboot <- apply(cover, 2, mean)
>  pboot.se <- sqrt(pboot * (1 - pboot) / nboot)
>  print(cbind(pboot, pboot.se))
     pboot  pboot.se
[1,]   0.8 0.1788854
[2,]   0.8 0.1788854
[3,]   1.0 0.0000000
[4,]   0.8 0.1788854
[5,]   1.0 0.0000000
[6,]   1.0 0.0000000
[7,]   1.0 0.0000000
> 
> ###################################################
> ### chunk number 38: doit (one mo time)
> ###################################################
> 
>  set.seed(42)
>  nboot <- 5
>  cover <- matrix(0, nboot, length(fit.hat))
>  for (iboot in 1:nboot) {
+      xstar <- raster(theta.hat, pred.orig, fam.orig, root.orig)
+      aout4star <- aster(xstar, root.orig, pred.orig, fam.orig,
+          modmat.orig, beta.hat, method = "nlm", check.analyticals = FALSE)
+      pout4star <- predict(aout4star, x.pred, root.pred, modmat.pred,
+          amat.orig, se.fit = TRUE)
+      upper <- pout4star$fit + crit.orig * pout4star$se.fit
+      lower <- pout4star$fit - crit.orig * pout4star$se.fit
+      cover[iboot, ] <- as.numeric(lower <= fit.hat & fit.hat <= upper)
+  }
>  pboot <- apply(cover, 2, mean)
>  pboot.se <- sqrt(pboot * (1 - pboot) / nboot)
>  print(cbind(pboot, pboot.se))
     pboot  pboot.se
[1,]   0.8 0.1788854
[2,]   0.8 0.1788854
[3,]   1.0 0.0000000
[4,]   0.8 0.1788854
[5,]   1.0 0.0000000
[6,]   1.0 0.0000000
[7,]   1.0 0.0000000
> 
> 
>