2 function(formula, data=NULL, origin=TRUE, nperm=999, method=NULL, silent=FALSE)
4 # This program computes a multiple linear regression and performs tests
5 # of significance of the equation parameters using permutations.
7 # origin=TRUE: the regression line can be forced through the origin. Testing
8 # the significance in that case requires a special permutation procedure.
10 # Permutation methods: raw data or residuals of full model
11 # Default method in regression through the origin: raw data
12 # Default method in ordinary multiple regression: residuals of full model
13 # - In ordinary multiple regression when m = 1: raw data
15 # Pierre Legendre, March 2009
17 if(!is.null(method)) method <- match.arg(method, c("raw", "residuals"))
18 if(is.null(method) & origin==TRUE) method <- "raw"
19 if(is.null(method) & origin==FALSE) method <- "residuals"
20 if(nperm < 0) stop("Incorrect value for 'nperm'")
22 ## From the formula, find the variables and the number of observations 'n'
23 toto <- lm(formula, data)
24 mf <- match.call(expand.dots = FALSE)
25 m <- match(c("formula", "data", "offset"), names(mf), 0)
27 mf$drop.unused.levels <- TRUE
28 mf[[1]] <- as.name("model.frame")
29 mf <- eval(mf, parent.frame())
30 var.names = colnames(mf) # Noms des variables
31 y <- as.matrix(mf[,1])
32 colnames(y) <- var.names[1]
33 X <- as.matrix(mf[,-1])
38 mm<- m # No. regression coefficients, possibly including the intercept
39 if(m == 1) method <- "raw"
40 if(nrow(X) != n) stop("Unequal number of rows in y and X")
43 if(!silent) cat("Regression through the origin",'\n')
46 if(!silent) cat("Multiple regression with estimation of intercept",'\n')
54 cat("Permutation method =",method,"data",'\n')
56 cat("Permutation method =",method,"of full model",'\n')
60 t.vec <- summary(reg)$coefficients[,3]
61 p.param.t <- summary(reg)$coefficients[,4]
62 df1 <- summary(reg)$fstatistic[[2]]
63 df2 <- summary(reg)$fstatistic[[3]]
64 F <- summary(reg)$fstatistic[[1]]
65 y.res <- summary(reg)$residuals
66 # b.vec <- summary(reg)$coefficients[,1]
67 # r.sq <- summary(reg)$r.squared
68 # adj.r.sq <- summary(reg)$adj.r.squared
69 # p.param.F <- pf(F, df1, df2, lower.tail=FALSE)
71 if(df1 < m) stop("\nCollinearity among the X variables. Check using 'lm'")
81 for(i in 1:nperm) # Permute raw data. Always use this method for F-test
83 if(origin) { # Regression through the origin
84 dia.bin <- diag((rbinom(n,1,0.5)*2)-1)
85 y.perm <- dia.bin %*% sample(y)
86 reg.perm <- lm(y.perm ~ 0 + X)
87 } else { # Multiple linear regression
89 reg.perm <- lm(y.perm ~ X)
92 # Permutation test of the F-statistic
93 F.perm <- summary(reg.perm)$fstatistic[1]
94 if(F.perm >= F) nGT.F <- nGT.F+1
96 # Permutation tests of the t-statistics: permute raw data
98 t.perm <- summary(reg.perm)$coefficients[,3]
99 if(nperm <= 5) cat(t.perm,'\n')
101 # One-tailed test in direction of sign
102 if(t.perm[j]*sign.t[j] >= t.vec[j]*sign.t[j]) nGT1.t[j] <- nGT1.t[j]+1
104 if( abs(t.perm[j]) >= abs(t.vec[j]) ) nGT2.t[j] <- nGT2.t[j]+1
109 if(method == "residuals") {
110 # Permute residuals of full model
112 if(origin) { # Regression through the origin
113 dia.bin <- diag((rbinom(n,1,0.5)*2)-1)
114 y.perm <- dia.bin %*% sample(y.res)
115 reg.perm <- lm(y.perm ~ 0 + X)
116 } else { # Multiple linear regression
117 y.perm <- sample(y.res,n)
118 reg.perm <- lm(y.perm ~ X)
121 # Permutation tests of the t-statistics: permute residuals
122 t.perm <- summary(reg.perm)$coefficients[,3]
123 if(nperm <= 5) cat(t.perm,'\n')
125 # One-tailed test in direction of sign
126 if(t.perm[j]*sign.t[j] >= t.vec[j]*sign.t[j]) nGT1.t[j] <- nGT1.t[j]+1
128 if( abs(t.perm[j]) >= abs(t.vec[j]) ) nGT2.t[j] <- nGT2.t[j]+1
132 # Compute the permutational probabilities
133 p.perm.F <- nGT.F/(nperm+1)
134 p.perm.t1 <- nGT1.t/(nperm+1)
135 p.perm.t2 <- nGT2.t/(nperm+1)
137 ### Do not test intercept by permutation of residuals in multiple regression
138 if(!origin & method=="residuals") {
139 if(silent) { # Note: silent==TRUE in simulation programs
140 p.perm.t1[1] <- p.perm.t2[1] <- 1
142 p.perm.t1[1] <- p.perm.t2[1] <- NA
148 a[3] <- sprintf("%2f",a[3])
149 if(!silent) cat("Computation time =",a[3]," sec",'\n')
153 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())
157 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())
161 class(out) <- "lmorigin"