1bgtest <- function(formula, order = 1, order.by = NULL, type = c("Chisq", "F"), data = list(), fill = 0) 2{ 3 dname <- paste(deparse(substitute(formula))) 4 5 if(!inherits(formula, "formula")) { 6 X <- if(is.matrix(formula$x)) 7 formula$x 8 else model.matrix(terms(formula), model.frame(formula)) 9 y <- if(is.vector(formula$y)) 10 formula$y 11 else model.response(model.frame(formula)) 12 } else { 13 mf <- model.frame(formula, data = data) 14 y <- model.response(mf) 15 X <- model.matrix(formula, data = data) 16 } 17 18 if(!is.null(order.by)) 19 { 20 if(inherits(order.by, "formula")) { 21 z <- model.matrix(order.by, data = data) 22 z <- as.vector(z[,ncol(z)]) 23 } else { 24 z <- order.by 25 } 26 X <- as.matrix(X[order(z),]) 27 y <- y[order(z)] 28 } 29 30 n <- nrow(X) 31 k <- ncol(X) 32 order <- 1:order 33 m <- length(order) 34 resi <- lm.fit(X,y)$residuals 35 36 Z <- sapply(order, function(x) c(rep(fill, length.out = x), resi[1:(n-x)])) 37 if(any(na <- !complete.cases(Z))) { 38 X <- X[!na, , drop = FALSE] 39 Z <- Z[!na, , drop = FALSE] 40 y <- y[!na] 41 resi <- resi[!na] 42 n <- nrow(X) 43 } 44 auxfit <- lm.fit(cbind(X,Z), resi) 45 46 cf <- auxfit$coefficients 47 vc <- chol2inv(auxfit$qr$qr) * sum(auxfit$residuals^2) / auxfit$df.residual 48 names(cf) <- colnames(vc) <- rownames(vc) <- c(colnames(X), paste("lag(resid)", order, sep = "_")) 49 50 switch(match.arg(type), 51 52 "Chisq" = { 53 bg <- n * sum(auxfit$fitted.values^2)/sum(resi^2) 54 p.val <- pchisq(bg, m, lower.tail = FALSE) 55 df <- m 56 names(df) <- "df" 57 }, 58 59 "F" = { 60 uresi <- auxfit$residuals 61 bg <- ((sum(resi^2) - sum(uresi^2))/m) / (sum(uresi^2) / (n-k-m)) 62 df <- c(m, n-k-m) 63 names(df) <- c("df1", "df2") 64 p.val <- pf(bg, df1 = df[1], df2 = df[2], lower.tail = FALSE) 65 }) 66 67 names(bg) <- "LM test" 68 RVAL <- list(statistic = bg, parameter = df, 69 method = paste("Breusch-Godfrey test for serial correlation of order up to", max(order)), 70 p.value = p.val, 71 data.name = dname, 72 coefficients = cf, 73 vcov = vc) 74 class(RVAL) <- c("bgtest", "htest") 75 return(RVAL) 76} 77 78vcov.bgtest <- function(object, ...) object$vcov 79df.residual.bgtest <- function(object, ...) if(length(df <- object$parameter) > 1L) df[2] else NULL 80 81