## see Greene (2003), Section 8.3.4, p.155
coxtest <- function(formula1, formula2, data = list())
{
## merge two models (if possible) and compute
## response y and regressor matrices X and Z
## 1. If formulas are available: build big model first
if(inherits(formula1, "formula") && inherits(formula2, "formula")) {
formula <- formula2
if(length(formula) > 2) formula[[2]] <- NULL
formula[[2]] <- as.call(list(as.name("+"), as.name("."), formula[[2]]))
formula <- update(formula1, formula)
mf <- model.frame(formula, data = data)
y <- model.response(mf)
X <- model.matrix(terms(formula1), mf)
Z <- model.matrix(terms(formula2), mf)
m1 <- deparse(formula1)
m2 <- deparse(formula2)
} else {
## 2. if not, go via row names
if(!inherits(formula1, "formula")) {
mf <- model.frame(formula1)
X <- if(is.matrix(formula1$x)) formula1$x else model.matrix(terms(formula1), mf)
y <- if(is.vector(formula1$y)) formula1$y else model.response(mf)
m1 <- deparse(formula(formula1))
} else {
mf <- model.frame(formula1, data = data)
y <- model.response(mf)
X <- model.matrix(formula1, data = data)
m1 <- deparse(formula1)
}
if(!inherits(formula2, "formula")) {
mf <- model.frame(formula2)
Z <- if(is.matrix(formula2$x)) formula2$x else model.matrix(terms(formula2), mf)
y2 <- if(is.vector(formula2$y)) formula2$y else model.response(mf)
m2 <- deparse(formula(formula2))
} else {
mf <- model.frame(formula2, data = data)
y2 <- model.response(mf)
Z <- model.matrix(formula2, data = data)
m2 <- deparse(formula2)
}
if(!(all(c(row.names(X) %in% row.names(Z), row.names(Z) %in% row.names(X))))) {
warning("models fitted on different subsets")
allnames <- row.names(X)[row.names(X) %in% row.names(Z)]
X <- X[allnames,]
Z <- Z[allnames,]
y <- y[allnames]
y2 <- y2[allnames]
}
if(!identical(y, y2)) warning("different dependent variables specified")
}
## for pretty printing
m1 <- paste(m1, collapse = "\n")
m2 <- paste(m2, collapse = "\n")
## Steps in Greene (2003) p.156-7
# 1.
x.fit <- lm.fit(X, y)
s2x <- mean(x.fit$residuals^2)
# 2.
z.fit <- lm.fit(Z, y)
s2z <- mean(z.fit$residuals^2)
# 3. / 5.
zx.fit <- lm.fit(Z, x.fit$fitted.values)
s2zx <- s2x + mean(zx.fit$residuals^2)
xz.fit <- lm.fit(X, z.fit$fitted.values)
s2xz <- s2z + mean(xz.fit$residuals^2)
# 4.
xzx.fit <- lm.fit(X, zx.fit$residuals)
s2xzx <- mean(xzx.fit$residuals^2)
zxz.fit <- lm.fit(Z, xz.fit$residuals)
s2zxz <- mean(zxz.fit$residuals^2)
# 6.
n <- nrow(X)
c01 <- n/2 * log(s2z/s2zx)
v01 <- n * (s2x/(s2zx^2)) * s2xzx
c10 <- n/2 * log(s2x/s2xz)
v10 <- n * (s2z/(s2xz^2)) * s2zxz
## Cox statistic for nonnested models
stat01 <- c01/sqrt(v01)
stat10 <- c10/sqrt(v10)
rval <- cbind(c(c01, c10), sqrt(c(v01, v10)), c(stat01, stat10),
2*pnorm(abs(c(stat01, stat10)), lower.tail=FALSE))
colnames(rval) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
rownames(rval) <- c("fitted(M1) ~ M2", "fitted(M2) ~ M1")
## put results together
title <- "Cox test\n"
topnote <- paste("Model ", 1:2,": ", c(m1, m2), sep="", collapse="\n")
rval <- structure(as.data.frame(rval), heading = c(title, topnote),
class = c("anova", "data.frame"))
return(rval)
}