Statistics 8931 (Geyer, Fall 2013) Example 2

R statements library(aster) data(echin2) vars <- unique(as.character(echin2$varb)) vars pred <- c(0, 1, 2, 3, 4, 5, 6, 6, 8, 8, 10) fam <- c(1, 1, 1, 1, 1, 1, 3, 1, 3, 1, 3) indoors <- grepl("lds", as.character(echin2$varb)) indoors <- as.numeric(indoors) outdoors <- 1 - indoors fit <- grepl("roct", as.character(echin2$varb)) fit <- as.numeric(fit) echin2 <- data.frame(echin2, indoors = indoors, outdoors = outdoors, fit = fit) layer <- gsub("[0-9]", "", as.character(echin2$varb)) layer <- as.factor(layer) echin2 <- data.frame(echin2, layer = layer) m <- model.matrix(resp ~ varb + layer : yearcross + indoors : flat + outdoors : layer : (row + posi) + fit : crosstype, data = echin2) outies <- grepl("row0|layerlds:outdoors|row10", colnames(m)) colnames(m)[outies] dim(m) m <- m[ , ! outies] dim(m) aout <- aster(resp ~ 0 + m, pred, fam, varb, id, root, data = echin2) summary(aout)

R statements library(aster) data(echin2) vars <- unique(as.character(echin2$varb)) vars pred <- c(0, 1, 2, 3, 4, 5, 6, 6, 8, 8, 10) fam <- c(1, 1, 1, 1, 1, 1, 3, 1, 3, 1, 3) indoors <- grepl("lds", as.character(echin2$varb)) indoors <- as.numeric(indoors) outdoors <- 1 - indoors fit <- grepl("roct", as.character(echin2$varb)) fit <- as.numeric(fit) echin2 <- data.frame(echin2, indoors = indoors, outdoors = outdoors, fit = fit) layer <- gsub("[0-9]", "", as.character(echin2$varb)) layer <- as.factor(layer) echin2 <- data.frame(echin2, layer = layer) m <- model.matrix(resp ~ varb + layer : yearcross + indoors : flat + outdoors : layer : (row + posi) + fit : crosstype, data = echin2) outies <- grepl("row0|layerlds:outdoors|row10", colnames(m)) colnames(m)[outies] dim(m) m <- m[ , ! outies] dim(m) aout <- aster(resp ~ 0 + m, pred, fam, varb, id, root, data = echin2) ########## above copied from box 1 ########## outies <- grepl("yearcross", colnames(m)) colnames(m)[outies] m2 <- m[ , ! outies] dim(m) dim(m2) aout2 <- aster(resp ~ 0 + m2, pred, fam, varb, id, root, data = echin2) outies <- grepl("crosstype", colnames(m2)) colnames(m2)[outies] m3 <- m2[ , ! outies] dim(m2) dim(m3) aout3 <- aster(resp ~ 0 + m3, pred, fam, varb, id, root, data = echin2) anova(aout3, aout2, aout) summary(aout2)

R statements library(aster) data(echin2) vars <- unique(as.character(echin2$varb)) vars pred <- c(0, 1, 2, 3, 4, 5, 6, 6, 8, 8, 10) fam <- c(1, 1, 1, 1, 1, 1, 3, 1, 3, 1, 3) indoors <- grepl("lds", as.character(echin2$varb)) indoors <- as.numeric(indoors) outdoors <- 1 - indoors fit <- grepl("roct", as.character(echin2$varb)) fit <- as.numeric(fit) echin2 <- data.frame(echin2, indoors = indoors, outdoors = outdoors, fit = fit) layer <- gsub("[0-9]", "", as.character(echin2$varb)) layer <- as.factor(layer) echin2 <- data.frame(echin2, layer = layer) m <- model.matrix(resp ~ varb + layer : yearcross + indoors : flat + outdoors : layer : (row + posi) + fit : crosstype, data = echin2) outies <- grepl("row0|layerlds:outdoors|row10", colnames(m)) colnames(m)[outies] dim(m) m <- m[ , ! outies] dim(m) aout <- aster(resp ~ 0 + m, pred, fam, varb, id, root, data = echin2) ########## above copied from box 1 ########## outies <- grepl("yearcross", colnames(m)) colnames(m)[outies] m2 <- m[ , ! outies] dim(m) dim(m2) aout2 <- aster(resp ~ 0 + m2, pred, fam, varb, id, root, data = echin2) outies <- grepl("crosstype", colnames(m2)) colnames(m2)[outies] m3 <- m2[ , ! outies] dim(m2) dim(m3) aout3 <- aster(resp ~ 0 + m3, pred, fam, varb, id, root, data = echin2) ########## above copied from box 2 ########## m4 <- model.matrix(resp ~ varb + layer : yearcross + indoors : flat + outdoors : layer : (row + posi) + fit : crosstype + indoors : crosstype, data = echin2) outies <- grepl("row0|layerlds:outdoors|row10", colnames(m4)) colnames(m4)[outies] dim(m4) m4 <- m4[ , ! outies] dim(m4) aout4 <- aster(resp ~ 0 + m4, pred, fam, varb, id, root, data = echin2) anova(aout3, aout2, aout4) outies <- grepl("fit:crosstype", colnames(m4)) colnames(m4)[outies] m5 <- m4[ , ! outies] dim(m4) dim(m5) aout5 <- aster(resp ~ 0 + m5, pred, fam, varb, id, root, data = echin2) anova(aout3, aout5, aout4)

R statements library(aster) data(echin2) vars <- unique(as.character(echin2$varb)) vars pred <- c(0, 1, 2, 3, 4, 5, 6, 6, 8, 8, 10) fam <- c(1, 1, 1, 1, 1, 1, 3, 1, 3, 1, 3) indoors <- grepl("lds", as.character(echin2$varb)) indoors <- as.numeric(indoors) outdoors <- 1 - indoors fit <- grepl("roct", as.character(echin2$varb)) fit <- as.numeric(fit) echin2 <- data.frame(echin2, indoors = indoors, outdoors = outdoors, fit = fit) layer <- gsub("[0-9]", "", as.character(echin2$varb)) layer <- as.factor(layer) echin2 <- data.frame(echin2, layer = layer) m <- model.matrix(resp ~ varb + layer : yearcross + indoors : flat + outdoors : layer : (row + posi) + fit : crosstype, data = echin2) outies <- grepl("row0|layerlds:outdoors|row10", colnames(m)) colnames(m)[outies] dim(m) m <- m[ , ! outies] dim(m) aout <- aster(resp ~ 0 + m, pred, fam, varb, id, root, data = echin2) ########## above copied from box 1 ########## outies <- grepl("yearcross", colnames(m)) colnames(m)[outies] m2 <- m[ , ! outies] dim(m) dim(m2) aout2 <- aster(resp ~ 0 + m2, pred, fam, varb, id, root, data = echin2) ########## above copied from box 2 ########## inies <- match(levels(echin2$crosstype), as.character(echin2$crosstype)) inies <- echin2$id[inies] inies newmodmat <- m2[echin2$id %in% inies, ] newmodmat[ , "indoors:flat1"] <- as.numeric(grepl("lds", rownames(newmodmat))) newmodmat[ , "indoors:flat2"] <- 0 newmodmat[ , "indoors:flat3"] <- 0 newmodmat[ , grep("row", colnames(newmodmat))] <- 0 newmodmat[ , grep("posi", colnames(newmodmat))] <- 0 nind <- nlevels(echin2$crosstype) nparm <- ncol(newmodmat) nnode <- length(newmodmat) / (nind * nparm) all.equal(nnode, as.integer(nnode)) newmodmat <- array(as.vector(newmodmat), c(nind, nnode, nparm)) newroot <- array(1, c(nind, nnode)) droppies <- gsub("m2", "", aout2$dropped) droppies outies <- match(droppies, colnames(m2)) outies newmodmat <- newmodmat[ , , - outies] dim(newmodmat) pout <- predict(aout2, root = newroot, modmat = newmodmat, se.fit = TRUE) pout.fit <- matrix(pout$fit, nrow = nind) pout.se <- matrix(pout$se.fit, nrow = nind) rownames(pout.fit) <- levels(echin2$crosstype) rownames(pout.se) <- levels(echin2$crosstype) colnames(pout.fit) <- vars colnames(pout.se) <- vars foo <- cbind(pout.fit[ , "lds3"], pout.se[ , "lds3"]) colnames(foo) <- c("fit", "se") foo