R/alr3.R

Defines functions anova.pod plot.pod plot.pod.lm predict.pod df.residual.pod fitted.pod formula.pod residuals.pod vcov.pod deviance.pod coef.pod summary.pod print.pod pod.formula pod.lm pod randomLinComb.lm randomLinComb.lm randomLinComb.default randomLinComb pureErrorAnova.lm pureErrorAnova

Documented in anova.pod coef.pod deviance.pod df.residual.pod fitted.pod formula.pod plot.pod plot.pod.lm pod pod.formula pod.lm predict.pod print.pod pureErrorAnova pureErrorAnova.lm randomLinComb randomLinComb.default randomLinComb.lm residuals.pod summary.pod vcov.pod

#  R version of the alr3 package  -- a few functions not moved to car

###########################################################################
# Chapter 5  pureErrorAnova
###########################################################################
# completely rewritten 2/24/2006
# renamed pureErrorAnova
pureErrorAnova <- function(mod){UseMethod("pureErrorAnova")}
pureErrorAnova.lm <- function(mod) {
 if (is.null(mod$model)) mod <- update(mod, model=TRUE)
 p <- dim(mod$model)[2] -1
 mod$model$Lack.of.Fit <- 
   factor(randomLinComb(model.matrix(mod), 101319853))
 aov1 <- anova(mod)
 #set.seed(save.seed) # restore random number seed
 if (length(levels(mod$model$Lack.of.Fit)) == length(mod$model$Lack.of.Fit)) 
  aov1 else {
  aov2 <- anova(lm(mod$model[ , 1]~mod$model$Lack.of.Fit, weights=weights(mod)))
  rrow <- dim(aov1)[1]
  aov2[1, 1] <- aov1[rrow, 1]-aov2[2, 1]
  aov2[1, 2] <- aov1[rrow, 2]-aov2[2, 2]
  aov2[1, 3] <- aov2[1, 2]/aov2[1, 1]
  aov1[1:(rrow-1), 4] <- aov1[1:(rrow-1), 3]/aov2[2, 3]
  aov2[1, 4] <- aov2[1, 3]/aov2[2, 3]
  row.names(aov2) <- c(" Lack of fit", " Pure Error")
  aov <- rbind(aov1, aov2)
  aov[ , 5] <- pf(aov[ , 4], aov[ , 1], aov2[2, 1], lower.tail=FALSE)
  aov
  }}

 
randomLinComb <- function(X, seed=NULL) {UseMethod("randomLinComb")}

randomLinComb.default <- function(X, seed=NULL) {
 if(!is.null(seed)) set.seed(seed)
 std <- function(x){ 
    s <- sd(x)
    if( s > 0) (x-mean(x))/s else x}
 as.vector(apply(X, 2, std)%*% as.vector(2*rnorm(dim(X)[2])-1) )
 }
 
randomLinComb.lm <- function(X, ...) {
  randomLinComb(model.matrix(X), ...)}
 
randomLinComb.lm <- function(X, seed=NULL) {
 if(is.null(X$model)) X <- update(X, model=TRUE)
 randomLinComb(X$model[ , -1], seed=seed)}


###########################################################################
# Chapter 6  Delta Method
###########################################################################
# in file deltaMethod.R in 'car'

##################################################################
# pod models; Cook and Weisberg (2004), American Statistician
##################################################################

# This code was written by Lexin Li.  It was modified and simplified 
# by S. Weisberg to work with standard R stuff in nls.

pod <- function(x, ...) {UseMethod("pod")}

pod.lm <- function(x, group, mean.function, control, ...) {
 call <- x$call
 call[[1]] <- as.name("pod")
 if(!missing(group)) call$group <- as.name(as.character(substitute(group)))
 if(!missing(mean.function)) call$mean.function <- 
     as.character(substitute(mean.function))
 if(!missing(control)) call$control <- 
     as.name(as.character(substitute(control)))
 eval(call)}

pod.formula <-
function (formula, data=sys.parent(), group, subset, weights=NULL, na.action, 
    mean.function=c("pod", "common", "parallel", "general"), 
    singular.ok = FALSE, contrasts = NULL, offset, control=nls.control(), ...) 
{ 
    op <- options(warn=-1)  # suppress warnings
    mod <- match.arg(mean.function)
    call <- match.call(expand.dots=FALSE)
    call$... <- NULL
    subset <- if (missing(subset)) NULL else call$subset
    g <- factor(eval(substitute(group), data))
    gname <- substitute(group)  
    gname <- as.character(if(length(gname) == 1) gname else gname[3])
    if (mod != "pod") { 
      if (mod == "common") {
       call$group <- call$mean.function <- NULL
       call[[1]] <- as.name("lm")
       ans <- eval(call, parent.frame())} else
     if (mod == "parallel") {
       assign(gname, g)
       call$formula <- update(as.formula(call$formula), 
                       as.formula(paste("~.+", gname, sep="")))
                        #as.formula("~.+g"))
       call$group <- call$mean.function <- NULL
       call[[1]] <- as.name("lm")
       ans <- eval(call)} else
     if (mod == "general") {
       assign(gname, g)
       call$formula <- update(as.formula(call$formula), 
                        as.formula(paste("~(.)*", gname, sep="")))
       call$group <- call$mean.function <- NULL
       call[[1]] <- as.name("lm")
       ans <- eval(call)}
     ans$pod.mean.function <- mod
     ans$group <- structure(if (!is.null(subset)) g[eval(subset)] else g,  
                            name=gname)
     class(ans) <- c("pod.lm", "lm")
     ans} else
    { # POD model
    l<-nlevels(g) 
    cl <- call
    cl$group <- cl$mean.function <- NULL
    cl[[1]] <- as.name("lm") 
    fit1 <- eval(cl)
    if(fit1$rank != length(fit1$coef))
     stop("Predictors are linearly related.  Try again after deleting predictors.")
    coef1 <- coef(fit1)
    p1 <- predict(fit1) - coef(fit1)[1]
# If subset is not NULL, find the subscripts of the cases used:
    if(!is.null(call$subset)){
       rows <- row.names(model.matrix(update(fit1, subset=NULL))) %in% 
               row.names(model.matrix(fit1))
       } else {
       rows <- rep(TRUE, dim(model.matrix(fit1))[1])}
       bign <- length(rows)
       ans <- rep(NA, bign)
       ans[rows] <- p1
       p1 <- ans
# update the formula in cl to fit parallel within group regression
    cl[[2]] <- update(as.formula(formula), as.formula("~ p1*g"))
# update cl$data
    if (is.null(cl$data)){
       cl$data <- data.frame(p1=p1, g=g)} else {
       cl$data <- data.frame(data, p1=p1, g=g)}
# fit parallel lines model  
    fit2 <- eval(cl)
    coef2 <- coef(fit2)
# group.data.subset excludes the subset
    group.data.subset <- 
              model.matrix(fit2$terms, fit2$model, contrasts)[ , 2+1:(l-1)]
    group.data.subset <- data.frame(group=group.data.subset)
# again, fix for non-null subset
    if(!is.null(subset)){
       group.data <- data.frame(matrix(NA, nrow=bign, 
          ncol=dim(group.data.subset)[2]))
       group.data[rows, ] <- group.data.subset} else {
       group.data <- group.data.subset}    
# construct new data with X, y, and G
    formula<-as.formula(formula)
    y.name<-as.character(attr(fit1$terms, "predvars")[2])
    x.name<-attr(fit1$terms, "term.labels")
    g.name<-paste(gname, levels(g)[2:l], sep="")
    names(group.data) <- g.name
    ncols <- length(x.name)+length(y.name)+length(g.name)
    p<-length(x.name)
# generate model
    para.name<-c("eta0", "eta1")
    form1<-paste("eta1 *", x.name[1])
    for(i in 2:p) {
       eta<-paste("eta", i, sep="")
       form1<-paste(form1, "+", eta, "*", x.name[i])
       para.name<-c(para.name, eta)  
    } 
   form2<-paste("eta0 +", form1)
   for(i in 1:(l-1)) {
      G<-g.name[i]
      th0<-paste("th0", i+1, sep="")
      th1<-paste("th1", i+1, sep="")
      form2<-paste(form2, "+", G, "* (", th0, "+", th1, "* (", form1, "))") 
      para.name<-c(para.name, c(th0, th1))
   } 
# New June 6, 2005
   if (is.null(weights))
      {form<-as.formula(paste(y.name, "~", form2))} else
      {wts <- substitute(weights)
       form<-as.formula(paste("sqrt(", wts, ")*", y.name, "~", 
             "sqrt(", wts, ")*(", form2, ")", sep=""))}
# End change
   start <- c(coef2[1], coef2[2]*coef1[-1], coef2[3:(l+1)], 
                coef2[(l+2):(2*l)]/coef2[2])
   names(start) <- c(paste("eta", 0:(length(x.name)), sep=""), 
                     paste("th0", 2:l, sep=""), paste("th1", 2:l, sep=""))
   if (is.null(call$data)){
        for (j in 1:length(g.name)) assign(g.name[j], group.data[ , j])
        obj.nl<-nls(form, start=start, subset=rows, 
            control=control, ...)} else {
        obj.nl<-nls(form, cbind(data, group.data), start=start, 
            subset=rows, control=control, ...)}  
# evaluate linear part, but only for the subset of cases used
   eta.name<-para.name[2:(p+1)]
   eta.value<-coef(obj.nl)[2:(p+1)]
   for(i in 1:p)
      assign(eta.name[i], eta.value[i])
   form1.expr<-parse(text=form1)
   envr <-
     if (is.null(call$data) & is.null(subset)){group.data} else{
      if(is.null(call$data)) group.data.subset else {
       if(is.null(subset)) as.data.frame(cbind(data, group.data)) else{
        as.data.frame(cbind(data[rows, ], group.data.subset))}}}
   if (!is.null(subset)) envr <- as.data.frame(envr[rows, ])
   linear.part<-eval(form1.expr, envr)
   ans <- NULL
   ans$nls.fit <- obj.nl
   ans$linear.part <- linear.part
   ans$group <- if (!is.null(subset)) g[rows] else g
   ans$call <- call
   class(ans)<-c("pod")
   options(op) # warnings ok again
   ans
}}

print.pod<-        function(x, ...){ print(x$nls.fit) }
summary.pod <-     function(object, ...){ summary(object$nls.fit, ...) }
coef.pod <-        function(object, ...){ coef(object$nls.fit, ...)}
deviance.pod <-    function(object, ...){ deviance(object$nls.fit, ...)} 
vcov.pod <-        function(object, ...){ vcov(object$nls.fit, ...)}
residuals.pod <-   function(object, ...){ residuals(object$nls.fit)}
formula.pod<-      function(x, ...){      formula(x$call)}
fitted.pod <-      function(object, ...){ fitted(object$nls.fit, ...)}
df.residual.pod <- function(object, ...){
                                 length(resid(object))-length(coef(object))}
predict.pod     <- function(object, ...){ predict(object$nls.fit)}

# Plot one dimensional models by group
plot.pod.lm <- function(x, colors=1:nlevels(x$group), 
      pch=1:nlevels(x$group), key="topleft", identify=FALSE, 
      xlab="Linear Predictor", ylab=as.character(c(formula(x)[[2]])), ...) {
  mean.function <- x$pod.mean.function
  if(mean.function == "general") stop("No 2D plot for the general pod model")
  g1 <- x$group
  g1.name <- attr(x$group, "name")
  levels(g1) <- 1:nlevels(x$group)
  g1 <- as.numeric(as.character(g1))
  gloc <- match("group", names(x$call))
  gname <- as.character(x$call[[gloc]])
# common regressions
  if(mean.function == "common"){
    plot(predict(x), x$model[ , 1], ylab=ylab, pch=pch[g1], col=colors[g1], 
    xlab=paste(xlab, ", ignore groups", sep=""), ...)
    abline(lm(x$model[ , 1]~predict(x)))}
# parallel regressions
  if(mean.function == "parallel"){
      c2 <- coef(x)
      xp <-0
      for (j in 2:(length(c2)-nlevels(x$group)+1)) 
         xp <- xp +c2[j] * model.matrix(x)[ , j]
      plot(xp, x$model[ , 1], pch=pch[g1], col=colors[g1], 
       ylab=paste(ylab, ", Groups = ", g1.name, sep=""), 
       xlab=paste(xlab,  ", parallel mean function", sep=""), ...)
      for (j in 1:nlevels(x$group))
       abline(if(j==1) c2[1] else c2[1]+c2[length(c2)-nlevels(x$group)+j], 1, 
           lty=j, col=colors[j])}
# key
  if (class(key) == "logical") {
   if (key == TRUE) {
      print("Click mouse on plot to locate the key, or press Escape")
      loc <-locator(n=1) 
      legend(loc[1], loc[2], legend = as.character(levels(x$group)), 
           lty=1:nlevels(x$group), col=colors[1:nlevels(x$group)], 
           pch=pch[1:nlevels(x$group)])}}
   else { 
      loc <- key
      legend(loc[1], loc[2], legend = as.character(levels(x$group)), 
           lty=1:nlevels(x$group), col=colors[1:nlevels(x$group)], 
           inset=0.01, pch=pch[1:nlevels(x$group)])}
# identify
  if(identify == TRUE){
      identify(xp, x$model[ , 1], row.names(x$model))}
  invisible()}

plot.pod <-function(x, colors=1:nlevels(x$group), 
  pch=1:nlevels(x$group), key="topleft", identify=FALSE, 
  xlab="Linear Predictor", ylab=as.character(c(formula(x)[[2]])), ...)
{ 
   yp<-residuals(x)+fitted(x)
   group <- x$group
   gloc <- match("group", names(x$call))
   xp<-x$linear.part
   g1 <- x$group
   levels(g1) <- 1:nlevels(group)
   g1 <- as.numeric(as.character(g1))
   gname <- as.character(x$call[[gloc]])
   plot(xp, yp, pch=pch[g1], col=colors[g1], 
    ylab=paste(ylab, ", Groups = ", as.character(x$call$group), sep=""), 
    xlab=paste(xlab,  ", pod mean function", sep=""), ...)
      for (j in 1:nlevels(group))
        {abline(lm(yp~xp, subset=g1==j), lty=j, col=colors[j])}
   if (class(key) == "logical") {
   if (key == TRUE) {
      print("Click mouse on plot to locate the key, or press Escape")
      loc <-locator(n=1) 
      legend(loc[1], loc[2], legend = as.character(levels(group)), 
           lty=1:nlevels(group), col=colors[1:nlevels(group)], 
           pch=pch[1:nlevels(group)])}} else
    { 
      loc <- key
      legend(loc[1], loc[2], legend = as.character(levels(group)), 
           lty=1:nlevels(group), col=colors[1:nlevels(group)], 
           inset=0.01, pch=pch[1:nlevels(group)])}
   invisible()
}

anova.pod <- 
function (object, scale = 0, test = "F", ...) 
{   
    m1 <- update(object, mean.function="common")
    m2 <- update(object, mean.function="parallel")
    m4 <- update(object, mean.function="general")
    objects <- list(m1, m2, object, m4)
    resdf  <- as.numeric(lapply(objects, df.residual))
    resdev <- as.numeric(lapply(objects, deviance))
    table <- data.frame(resdf, resdev, c(NA, -diff(resdf)), c(NA, 
        -diff(resdev)))
    variables <- c("1: common", "2: parallel", "3: pod", "4: pod + 2fi")
    dimnames(table) <- list(variables, c("Res.Df", "RSS", "Df", 
        "Sum of Sq"))
    title <- paste("POD Analysis of Variance Table for ", 
            deparse(formula(objects[[1]])[[2]]), ", grouped by ", 
            as.character(object$call$group), "\n" , sep="")
    topnote <- c(paste("1: ", deparse(formula(objects[[1]])), sep=""), 
                 paste("2: ", deparse(formula(objects[[2]])), sep=""), 
                 paste("3: ", deparse(formula(object$nls.fit)), sep=""), 
                 paste("4: ", deparse(formula(objects[[4]])), sep=""))
    if (!is.null(test)) {
        bigmodel <- order(resdf)[1]
        scale <- if (scale > 0) 
            scale
        else resdev[bigmodel]/resdf[bigmodel]
        table <- stat.anova(table = table, test = test, scale = scale, 
            df.scale = resdf[bigmodel], n = length(objects[bigmodel$residuals]))
    }
    structure(table, heading = c(title, topnote), class = c("anova", 
        "data.frame"))
}


#######################################################################
# Chapter 7, Transformations
#######################################################################
# Transformations --- all moved to 'car'

####################################################################
#inv.tran.plot now 'invTranPlot' in 'car'

##########################################################################
# replaced by 'powerTransform' in 'car'
    
    
################################################################################
#    Chapter 8
################################################################################
# Test for curvature in a residual plot (Chapter 8)
# Residual plots, and curvature tests in 'car'
  
#resid.curv.test = residCurvTest in 'car'

#resplot = residualPlot in 'car'

#############################################
# marginal model plots    Rev 10/30/07
# mmps = mmps in 'car'

################################################################
# Chapter 9
################################################################

#inf.index = infIndexPlot in 'car'

#outlier.t.test = outlierTest in 'car'

Try the alr3 package in your browser

Any scripts or data that you put into this service are public.

alr3 documentation built on Nov. 20, 2017, 3:02 a.m.