# S3 method for class "mlm"
influence.mlm <- function(model, m=1, do.coef=TRUE, ...)
	mlm.influence(model, m=m, do.coef = do.coef, ...)


# main computation a la lm.influence
mlm.influence <-
function (model, m=1, do.coef = TRUE, ...) 
{

	# helper functions
  vec <- function(M) {
  	R <- matrix(M, ncol=1)
  	if (is.vector(M)) return(R)
  	nn<-expand.grid(dimnames(M))[,2:1]
  	rownames(R) <- apply(as.matrix(nn), 1, paste, collapse=":")
  	R
  	}

	X <- model.matrix(model)
	data <- model.frame(model)
	Y <- as.matrix(model.response(data))
	r <- ncol(Y)
	n <- nrow(X)
	p <- ncol(X)
	labels <- rownames(X)
	call <- model$call
	
	B <- coef(model)
	E <- residuals(model)
	XPXI <- solve(crossprod(X))
	EPEI <- solve(crossprod(E))
	vB <- vec(t(B))
	S <- crossprod(E)/(n-p);
	V <- solve(p*S) %x% crossprod(X)
	
	subsets <- t(combn(n, m))
	nsub <- nrow(subsets) 
	

	Beta <- as.list(rep(0, nsub))
	R <- L <- H <- Q <- as.list(rep(0, nsub))
	CookD <- as.vector(rep(0, nsub))

	# FIXME: naive, direct computations can use qr() or update()
	for (i in seq(nsub)) {
		I <- c(subsets[i,])
		rows <- which(!(1:n) %in% I)
		XI <- X[rows,]
		YI <- Y[rows,]
		BI <- solve(crossprod(XI)) %*% t(XI) %*% YI
		EI <- (Y - X %*% BI)[I, , drop=FALSE]
		CookD[i] <- t(vec(B - BI)) %*% V %*% vec(B - BI)
		H[[i]] <- X[I, , drop=FALSE] %*% XPXI %*% t(X[I, , drop=FALSE])
		Q[[i]] <- EI %*% EPEI %*% t(EI)
		if (do.coef) Beta[[i]] <- BI
		L[[i]] <- if(m==1) H[[i]] / (1-H[[i]]) else H[[i]] %*% solve(diag(m) - H[[i]])
		R[[i]] <- if(m==1) Q[[i]] / (1-H[[i]]) 
					else {
					IH <- mpower(diag(m)-H[[i]], -1/2)
					IH %*% Q[[i]] %*% IH
					}
	}
	if(m==1) {
		H <- unlist(H)
		Q <- unlist(Q)
		L <- unlist(L)
		R <- unlist(R)
		subsets <- c(subsets)
	}
	result <- list(m=m, H=H, Q=Q, CookD=CookD, L=L, R=R, subsets=subsets, labels=labels, call=call)
	if (do.coef) result <- c(result, list(Beta=Beta))
	class(result) <-"inflmlm"
	result
} 

#####################
# extractor functions -- trying to find a way to be consistent with the .lm versions
hatvalues.mlm <- function(model, m=1, infl=mlm.influence(model, m=m, do.coef = FALSE), ...) 
{
    hat <- infl$H
    m <- infl$m
    names(hat) <- if(m==1) infl$subsets else apply(infl$subsets,1, paste, collapse=',')
    hat
}

cooks.distance.mlm <-
function (model, infl = mlm.influence(model, do.coef = FALSE), ...) 
{
    cookd <- infl$CookD
    m <- infl$m
    names(cookd) <- if(m==1) infl$subsets else apply(infl$subsets,1, paste, collapse=',')
    cookd
}

TESTME <- TRUE
if(TESTME) {
library(heplots)
Rohwer2 <- subset(Rohwer, subset=group==2)
rownames(Rohwer2)<- 1:nrow(Rohwer2)
Rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ n+s+ns+na+ss, data=Rohwer2)

## this doesn't work, as I would like it to, calling the hatvalues.mlm method, but passing m=2:
hatvalues(Rohwer.mod, m=2)
## Error in UseMethod("hatvalues") :
##  no applicable method for 'hatvalues' applied to an object of class "c('double', 'numeric')"

#I don't understand why this doesn't just call hatvalues.mlm, since Rohwer.mod is of class "mlm".

# These work -- calling hatvalues.mlm explicitly, or passing the infl= argument with the call to
# mlm.influence
hatvalues.mlm(Rohwer.mod, m=2)
hatvalues(Rohwer.mod, infl=mlm.influence(Rohwer.mod,m=2))

}