'lmorigin' <- 
function(formula, data=NULL, origin=TRUE, nperm=999, method=NULL, silent=FALSE)
#
# This program computes a multiple linear regression and performs tests 
# of significance of the equation parameters using permutations. 
# 
# origin=TRUE: the regression line can be forced through the origin. Testing 
# the significance in that case requires a special permutation procedure.
#
# Permutation methods: raw data or residuals of full model
#    Default method in regression through the origin: raw data
#    Default method in ordinary multiple regression: residuals of full model
#    - In ordinary multiple regression when m = 1: raw data
#
#         Pierre Legendre, March 2009
{
if(!is.null(method)) method <- match.arg(method, c("raw", "residuals"))
if(is.null(method) & origin==TRUE) method <- "raw"
if(is.null(method) & origin==FALSE) method <- "residuals"
if(nperm < 0) stop("Incorrect value for 'nperm'")
	
## From the formula, find the variables and the number of observations 'n'
toto <- lm(formula, data)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "offset"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
var.names = colnames(mf)   # Noms des variables
y <- as.matrix(mf[,1])
colnames(y) <- var.names[1]
X <- as.matrix(mf[,-1])
n <- nrow(mf)
m <- ncol(X)

a <- system.time({
mm<- m           # No. regression coefficients, possibly including the intercept
if(m == 1) method <- "raw"
if(nrow(X) != n) stop("Unequal number of rows in y and X")

if(origin) {
	if(!silent) cat("Regression through the origin",'\n')
	reg <- lm(y ~ 0 + X)
	} else {
	if(!silent) cat("Multiple regression with estimation of intercept",'\n')
	reg <- lm(y ~ X)	
	mm <- mm+1
	}

if(!silent) {
	if(nperm > 0) {
		if(method == "raw") {
			cat("Permutation method =",method,"data",'\n')
			} else {
			cat("Permutation method =",method,"of full model",'\n')
			}
		}
	}
t.vec      <- summary(reg)$coefficients[,3]
p.param.t  <- summary(reg)$coefficients[,4]
df1        <- summary(reg)$fstatistic[[2]]
df2        <- summary(reg)$fstatistic[[3]]
F          <- summary(reg)$fstatistic[[1]]
y.res      <- summary(reg)$residuals
# b.vec      <- summary(reg)$coefficients[,1]
# r.sq       <- summary(reg)$r.squared
# adj.r.sq   <- summary(reg)$adj.r.squared
# p.param.F  <- pf(F, df1, df2, lower.tail=FALSE)

if(df1 < m) stop("\nCollinearity among the X variables. Check using 'lm'")

# Permutation tests
if(nperm > 0) {

	nGT.F  <- 1
	nGT1.t <- rep(1,mm)
	nGT2.t <- rep(1,mm)
	sign.t <- sign(t.vec)

	for(i in 1:nperm)  # Permute raw data. Always use this method for F-test
		{
		if(origin) {   # Regression through the origin
			dia.bin <- diag((rbinom(n,1,0.5)*2)-1)
			y.perm <- dia.bin %*% sample(y)
			reg.perm <- lm(y.perm ~ 0 + X)
			} else {   # Multiple linear regression
			y.perm <- sample(y,n)
			reg.perm <- lm(y.perm ~ X)
			}
	
		# Permutation test of the F-statistic
		F.perm <- summary(reg.perm)$fstatistic[1]
		if(F.perm >= F) nGT.F <- nGT.F+1

		# Permutation tests of the t-statistics: permute raw data
		if(method == "raw") {
			t.perm <- summary(reg.perm)$coefficients[,3]
			if(nperm <= 5) cat(t.perm,'\n')
			for(j in 1:mm) {
				# One-tailed test in direction of sign
				if(t.perm[j]*sign.t[j] >= t.vec[j]*sign.t[j]) nGT1.t[j] <- nGT1.t[j]+1
				# Two-tailed test
				if( abs(t.perm[j]) >= abs(t.vec[j]) ) nGT2.t[j] <- nGT2.t[j]+1
				}
			}
		}
		
	if(method == "residuals") {
	# Permute residuals of full model
		for(i in 1:nperm) {
			if(origin) {   # Regression through the origin
				dia.bin <- diag((rbinom(n,1,0.5)*2)-1)
				y.perm <- dia.bin %*% sample(y.res)
				reg.perm <- lm(y.perm ~ 0 + X)
				} else {   # Multiple linear regression
				y.perm <- sample(y.res,n)
				reg.perm <- lm(y.perm ~ X)
				}
	
			# Permutation tests of the t-statistics: permute residuals
			t.perm <- summary(reg.perm)$coefficients[,3]
			if(nperm <= 5) cat(t.perm,'\n')
			for(j in 1:mm) {
				# One-tailed test in direction of sign
				if(t.perm[j]*sign.t[j] >= t.vec[j]*sign.t[j]) nGT1.t[j] <- nGT1.t[j]+1
				# Two-tailed test
				if( abs(t.perm[j]) >= abs(t.vec[j]) ) nGT2.t[j] <- nGT2.t[j]+1
				}
			}
		}
	# Compute the permutational probabilities
	p.perm.F <- nGT.F/(nperm+1)
	p.perm.t1 <- nGT1.t/(nperm+1)
	p.perm.t2 <- nGT2.t/(nperm+1)

	### Do not test intercept by permutation of residuals in multiple regression
	if(!origin & method=="residuals") {
		if(silent) {   # Note: silent==TRUE in simulation programs
			p.perm.t1[1] <- p.perm.t2[1] <- 1
			} else {
			p.perm.t1[1] <- p.perm.t2[1] <- NA
			}
		}
	}

})
a[3] <- sprintf("%2f",a[3])
if(!silent) cat("Computation time =",a[3]," sec",'\n')
#
if(nperm == 0) {

	out <- list(reg=reg, p.param.t.2tail=p.param.t, p.param.t.1tail=p.param.t/2, origin=origin, nperm=nperm, var.names=var.names, call=match.call())

	} else {

	out <- list(reg=reg, p.param.t.2tail=p.param.t, p.param.t.1tail=p.param.t/2, p.perm.t.2tail=p.perm.t2, p.perm.t.1tail=p.perm.t1, p.perm.F=p.perm.F, origin=origin, nperm=nperm, method=method, var.names=var.names, call=match.call())

	}
#
class(out) <- "lmorigin"
out
}