##########################
##### LOAD LIBRARIES #####
##########################

library (lattice)
library(Design)
library(Hmisc)

##################################
##### DESCRIPTIVE STATISTICS #####
##################################

# Standard error
se <- function (variable) {
	st.error = sd(na.omit(variable))/sqrt(length(na.omit(variable)))
	return(st.error)
}

# Z-scores
z.scores <- function(vector,mean=NULL,sd=NULL) {
	if (length(levels(as.factor(mean))) == 0) {
		mean=mean(vector)
	}
	if (length(levels(as.factor(sd))) == 0) {
		sd=sd(vector)
	}
	zscores = (vector-mean)/sd
	return(zscores)
}

###############################
##### LOGISTIC REGRESSION #####
###############################

odds <- function(p) {p/(1-p)}
logit <- function (x) {log(odds(x))}
invlogit <- function (x) {1/(1+exp(-x))}

#######################
##### EFFECT SIZE #####
#######################

##### R-SQUARED: Regression, and any General Linear Models (including ANOVAS)
r.squared <- function (model) {
	r.sq = cor(fitted(model),model$model[,1])^2
	return(r.sq)
}

##### ETA-SQUARED: ANOVA (eta-squared is a special case of R-squared)
eta.squared <- function (model.aov) {
	nvariables = dim(model.aov$model)[2]
	nobservations = dim(model.aov$model)[1]
	ncontrasts = 0
	for (i in 2:nvariables) {
		ncontrasts = ncontrasts + length(as.factor(levels(as.factor(model.aov$model[,i]))))-1
	}
	endeffects = ncontrasts+1
	startresiduals = ncontrasts+2
	sseffects = sum(model.aov$effects[2:endeffects]^2)
	ssresiduals = sum(model.aov$effects[startresiduals:nobservations]^2)
	eta.sq = sseffects / (sseffects + ssresiduals)
	return (eta.sq)
}

##### OMEGA-SQUARED: ANOVA (a corrected effect-size index)
omega.squared <- function (model) {
	if (class(model)[1]=="aov") {
		df.factors=0
		for (i in 2:dim(model$model)[2]) {
			df.factors=df.factors+length(levels(model$model[,i]))-1
		}
	}
	if (class(model)[1]=="lm") {
		df.factors=summary(model)$fstatistic[2]
	}
	
	nvariables = dim(model$model)[2]
	nobservations = dim(model$model)[1]
	ncontrasts = 0
	for (i in 2:nvariables) {
		ncontrasts = ncontrasts + length(as.factor(levels(as.factor(model$model[,i]))))-1
	}
	endeffects = ncontrasts+1
	startresiduals = ncontrasts+2
	sseffects = sum(model$effects[2:endeffects]^2)
	ssresiduals = sum(model$effects[startresiduals:nobservations]^2)
	
	omega.sq = (sseffects-(df.factors*(ssresiduals/model$df.residual)))/(sseffects+ssresiduals+(ssresiduals/model$df.residual))
	attr(omega.sq,"names") <- NULL
	return(omega.sq)
}

##### F-SQUARED: Regression. Ratio: Proportion of Variance accounted for by the effects, over Proportion of Variance not accounted for by the effects.
f.squared <- function (model) {
	f.sq = r.squared(model)/(1-r.squared(model))
	return(f.sq)
}

##### F: ANOVA. Standard deviation of standardized means. Just calculate the sqrt of f.squared

##### COHEN'S d: T-test.
cohen.d <- function (variable1, variable2) {
	cohens.d = (mean(variable1)-mean(variable2))/sqrt(((length(variable1)-1)*sd(variable1)^2+(length(variable2)-1)*sd(variable2)^2)/(length(variable1)+length(variable2)))
	return(cohens.d)
}

##### OMEGA: Chi-squared test.
omega <- function (matrix) {
	w = sqrt(chisq.test(matrix)$statistic/sum(matrix))
	attr(w,"names")<- NULL
	return (w)
}

##### PARTIAL R-SQUARED. Calculate the increase in the variance accounted for by the addition of successive factors in a model.
# In partial R2, the calculated ratio is: (R2.increase) over (1 - R2.of.previous.(subset).model).
# It is interpreted as the variance accounted for by the new factor when the previous factor(s) are constant or controled.
partial.r.squared <- function (model, method="all") {
	nfactors = length(names(model$model)) -1
	p.values = anova(model)[1:nfactors,5]
	p.order=order(p.values)
	nsignif001 = length(p.values[p.values <= 0.001])
	nsignif01 = length(p.values[p.values <= 0.01])
	nsignif05 = length(p.values[p.values <= 0.05])
	
	if (method == "all") {
		nused = nfactors
	}
	if (method == "significant") {
		nused = length(p.values[p.values <= 0.05])
	}
	
	stars = c(rep(" (***)",nsignif001),rep(" (**)",nsignif01-nsignif001),rep(" (*)",nsignif05-nsignif01),rep("",nused-nsignif05))

	used.factors = p.order[1:nused]
	
	endfactors = nfactors+1
	factor.names = names(model$model)[2:endfactors]
	factor.names = factor.names[used.factors]
	factor.names = paste (factor.names,stars,sep="")
	used.factors = used.factors +1

	partial.r2 = NULL

	if (nused >= 1) {
		lm.1 = lm (model$model[,1] ~ model$model[,used.factors[1]])
		part.r2 = r.squared (lm.1)
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 2) {
		lm.2 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]])
		part.r2 = (r.squared(lm.2) - r.squared(lm.1)) / (1 - r.squared(lm.1))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 3) {
		lm.3 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]])
		part.r2 = (r.squared(lm.3) - r.squared(lm.2)) / (1 - r.squared(lm.2))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 4) {
		lm.4 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]])
		part.r2 = (r.squared(lm.4) - r.squared(lm.3)) / (1 - r.squared(lm.3))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 5) {
		lm.5 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]])
		part.r2 = (r.squared(lm.5) - r.squared(lm.4)) / (1 - r.squared(lm.4))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 6) {
		lm.6 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]])
		part.r2 = (r.squared(lm.6) - r.squared(lm.5)) / (1 - r.squared(lm.5))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 7) {
		lm.7 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]])
		part.r2 = (r.squared(lm.7) - r.squared(lm.6)) / (1 - r.squared(lm.6))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 8) {
		lm.8 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]])
		part.r2 = (r.squared(lm.8) - r.squared(lm.7)) / (1 - r.squared(lm.7))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 9) {
		lm.9 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]])
		part.r2 = (r.squared(lm.9) - r.squared(lm.8)) / (1 - r.squared(lm.8))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 10) {
		lm.10 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]])
		part.r2 = (r.squared(lm.10) - r.squared(lm.9)) / (1 - r.squared(lm.9))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 11) {
		lm.11 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]])
		part.r2 = (r.squared(lm.11) - r.squared(lm.10)) / (1 - r.squared(lm.10))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 12) {
		lm.12 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]])
		part.r2 = (r.squared(lm.12) - r.squared(lm.11)) / (1 - r.squared(lm.11))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 13) {
		lm.13 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]])
		part.r2 = (r.squared(lm.13) - r.squared(lm.12)) / (1 - r.squared(lm.12))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 14) {
		lm.14 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]] + model$model[,used.factors[14]])
		part.r2 = (r.squared(lm.14) - r.squared(lm.13)) / (1 - r.squared(lm.13))
		partial.r2 = c(partial.r2, part.r2)
	}

	if (nused >= 15) {
		lm.15 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]] + model$model[,used.factors[14]] + model$model[,used.factors[15]])
		part.r2 = (r.squared(lm.15) - r.squared(lm.14)) / (1 - r.squared(lm.14))
		partial.r2 = c(partial.r2, part.r2)
	}

	partial.r2 = round(partial.r2,digits=4)
	attr(partial.r2,"names") <- factor.names
	return (partial.r2)
}

##### PARTIAL ETA-SQUARED. Special case of partial R2 for ANOVA.
partial.eta.squared <- function (model,method="all") {
	if (method == "all") {
		partial.e2 = partial.r.squared (model)
	}
	if (method == "significant") {
		partial.e2 = partial.r.squared (model,method="significant")
	}
	return (partial.e2)
}

##### SEMIPARTIAL R-SQUARED. Calculate the increase in the variance accounted for by the addition of successive factors in a model.
# In partial R2, the calculated ratio is: (R2.increase) over (1 - R2.of.total.(current).model).
# It is interpreted as the variance accounted for by the new factor when the previous factor(s) is partialled out from the new one, but not from the dependent variable.
semipartial.r.squared <- function (model, method="all") {
	nfactors = length(names(model$model)) -1
	p.values = anova(model)[1:nfactors,5]
	p.order=order(p.values)
	nsignif001 = length(p.values[p.values <= 0.001])
	nsignif01 = length(p.values[p.values <= 0.01])
	nsignif05 = length(p.values[p.values <= 0.05])
	
	if (method == "all") {
		nused = nfactors
	}
	if (method == "significant") {
		nused = length(p.values[p.values <= 0.05])
	}
	
	stars = c(rep(" (***)",nsignif001),rep(" (**)",nsignif01-nsignif001),rep(" (*)",nsignif05-nsignif01),rep("",nused-nsignif05))

	used.factors = p.order[1:nused]
	
	endfactors = nfactors+1
	factor.names = names(model$model)[2:endfactors]
	factor.names = factor.names[used.factors]
	factor.names = paste (factor.names,stars,sep="")
	used.factors = used.factors +1

	semipartial.r2 = NULL

	if (nused >= 1) {
		lm.1 = lm (model$model[,1] ~ model$model[,used.factors[1]])
		semipart.r2 = r.squared (lm.1)
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 2) {
		lm.2 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]])
		semipart.r2 = (r.squared(lm.2) - r.squared(lm.1)) / (1 - r.squared(lm.2))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 3) {
		lm.3 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]])
		semipart.r2 = (r.squared(lm.3) - r.squared(lm.2)) / (1 - r.squared(lm.3))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 4) {
		lm.4 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]])
		semipart.r2 = (r.squared(lm.4) - r.squared(lm.3)) / (1 - r.squared(lm.4))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 5) {
		lm.5 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]])
		semipart.r2 = (r.squared(lm.5) - r.squared(lm.4)) / (1 - r.squared(lm.5))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 6) {
		lm.6 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]])
		semipart.r2 = (r.squared(lm.6) - r.squared(lm.5)) / (1 - r.squared(lm.6))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 7) {
		lm.7 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]])
		semipart.r2 = (r.squared(lm.7) - r.squared(lm.6)) / (1 - r.squared(lm.7))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 8) {
		lm.8 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]])
		semipart.r2 = (r.squared(lm.8) - r.squared(lm.7)) / (1 - r.squared(lm.8))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 9) {
		lm.9 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]])
		semipart.r2 = (r.squared(lm.9) - r.squared(lm.8)) / (1 - r.squared(lm.9))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 10) {
		lm.10 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]])
		semipart.r2 = (r.squared(lm.10) - r.squared(lm.9)) / (1 - r.squared(lm.10))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 11) {
		lm.11 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]])
		semipart.r2 = (r.squared(lm.11) - r.squared(lm.10)) / (1 - r.squared(lm.11))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 12) {
		lm.12 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]])
		semipart.r2 = (r.squared(lm.12) - r.squared(lm.11)) / (1 - r.squared(lm.12))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 13) {
		lm.13 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]])
		semipart.r2 = (r.squared(lm.13) - r.squared(lm.12)) / (1 - r.squared(lm.13))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 14) {
		lm.14 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]] + model$model[,used.factors[14]])
		semipart.r2 = (r.squared(lm.14) - r.squared(lm.13)) / (1 - r.squared(lm.14))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	if (nused >= 15) {
		lm.15 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]] + model$model[,used.factors[14]] + model$model[,used.factors[15]])
		semipart.r2 = (r.squared(lm.15) - r.squared(lm.14)) / (1 - r.squared(lm.15))
		semipartial.r2 = c(semipartial.r2, semipart.r2)
	}

	semipartial.r2 = round(semipartial.r2,digits=4)
	attr(semipartial.r2,"names") <- factor.names
	return (semipartial.r2)
}

##### ACCUMULATED R-SQUARED. Calculate R2 for successive nested models.
# It doesn't calculate the increase, but the final R2 of the total model.
accumulated.r.squared <- function (model, method="all") {
	nfactors = length(names(model$model)) -1
	p.values = anova(model)[1:nfactors,5]
	p.order=order(p.values)
	nsignif001 = length(p.values[p.values <= 0.001])
	nsignif01 = length(p.values[p.values <= 0.01])
	nsignif05 = length(p.values[p.values <= 0.05])
	
	if (method == "all") {
		nused = nfactors
	}
	if (method == "significant") {
		nused = length(p.values[p.values <= 0.05])
	}
	
	stars = c(rep(" (***)",nsignif001),rep(" (**)",nsignif01-nsignif001),rep(" (*)",nsignif05-nsignif01),rep("",nused-nsignif05))

	used.factors = p.order[1:nused]
	
	endfactors = nfactors+1
	factor.names = names(model$model)[2:endfactors]
	factor.names = factor.names[used.factors]
	factor.names = paste (factor.names,stars,sep="")
	used.factors = used.factors +1

	accum.r2 = NULL

	if (nused >= 1) {
		lm.1 = lm (model$model[,1] ~ model$model[,used.factors[1]])
		accum.r2 = c(accum.r2, r.squared(lm.1))
	}

	if (nused >= 2) {
		lm.2 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]])
		accum.r2 = c(accum.r2, r.squared(lm.2))
	}

	if (nused >= 3) {
		lm.3 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]])
		accum.r2 = c(accum.r2, r.squared(lm.3))
	}

	if (nused >= 4) {
		lm.4 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]])
		accum.r2 = c(accum.r2, r.squared(lm.4))
	}

	if (nused >= 5) {
		lm.5 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]])
		accum.r2 = c(accum.r2, r.squared(lm.5))
	}

	if (nused >= 6) {
		lm.6 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]])
		accum.r2 = c(accum.r2, r.squared(lm.6))
	}

	if (nused >= 7) {
		lm.7 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]])
		accum.r2 = c(accum.r2, r.squared(lm.7))
	}

	if (nused >= 8) {
		lm.8 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]])
		accum.r2 = c(accum.r2, r.squared(lm.8))
	}

	if (nused >= 9) {
		lm.9 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]])
		accum.r2 = c(accum.r2, r.squared(lm.9))
	}

	if (nused >= 10) {
		lm.10 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]])
		accum.r2 = c(accum.r2, r.squared(lm.10))
	}

	if (nused >= 11) {
		lm.11 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]])
		accum.r2 = c(accum.r2, r.squared(lm.11))
	}

	if (nused >= 12) {
		lm.12 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]])
		accum.r2 = c(accum.r2, r.squared(lm.12))
	}

	if (nused >= 13) {
		lm.13 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]])
		accum.r2 = c(accum.r2, r.squared(lm.13))
	}

	if (nused >= 14) {
		lm.14 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]] + model$model[,used.factors[14]])
		accum.r2 = c(accum.r2, r.squared(lm.14))
	}

	if (nused >= 15) {
		lm.15 = lm (model$model[,1] ~ model$model[,used.factors[1]] + model$model[,used.factors[2]] + model$model[,used.factors[3]] + model$model[,used.factors[4]] + model$model[,used.factors[5]] + model$model[,used.factors[6]] + model$model[,used.factors[7]] + model$model[,used.factors[8]] + model$model[,used.factors[9]] + model$model[,used.factors[10]] + model$model[,used.factors[11]] + model$model[,used.factors[12]] + model$model[,used.factors[13]] + model$model[,used.factors[14]] + model$model[,used.factors[15]])
		accum.r2 = c(accum.r2, r.squared(lm.15))
	}

	accum.r2 = round(accum.r2,digits=4)
	attr(accum.r2,"names") <- factor.names
	return (accum.r2)
}

####################
##### GRAPHICS #####
####################

##### BAR (COLUMN) PLOT WITH STANDARD-ERROR BARS
bar.and.error <- function (variable1, variable2, variable3=NULL, xlim=NULL,legend.text=NULL,args.legend=NULL,legend.margin=1,col=NULL,main=NULL,sub=NULL,xlab=NULL,ylab=NULL) {
	levels2=as.factor(unique(as.factor(variable2)))
	nlevels2=length(levels2)
	levels3=as.factor(unique(as.factor(variable3)))
	nlevels3=length(levels3)
	mean.matrix=NULL
	se.matrix=NULL
	abcissa.matrix=NULL
	if (nlevels3 == 0) {
		current.abcissa = 0.7
		for (i in 1:nlevels2) {
			mean.i = mean (variable1 [as.factor(variable2)==levels2[i]])
			mean.matrix = c(mean.matrix,mean.i)
			se.i = se (variable1 [as.factor(variable2)==levels2[i]])
			se.matrix = c(se.matrix,se.i)
			abcissa.matrix = c(abcissa.matrix,current.abcissa)
			current.abcissa = current.abcissa + 1.2
		}
	}
	else {
		legend.text=as.character(unique(as.factor(variable3)))
		xlim=c(0,0.5+nlevels2*(nlevels3+1)+legend.margin)
		for (i in 1:nlevels2) {
			mean.vector=NULL
			se.vector=NULL
			abcissa.vector=NULL
			current.abcissa=1.5+((nlevels3+1)*(i-1))
			for (j in 1:nlevels3) {
				current.abcissa=current.abcissa+j-1
				mean.ij = mean (variable1 [as.factor(variable2)==levels2[i] & as.factor(variable3)==levels3[j]])
				mean.vector = c(mean.vector,mean.ij)
				se.ij = se (variable1 [as.factor(variable2)==levels2[i] & as.factor(variable3)==levels3[j]])
				se.vector = c(se.vector,se.ij)
				abcissa.vector = c(abcissa.vector,current.abcissa)
			}
			mean.matrix=cbind(mean.matrix,mean.vector)
			se.matrix = cbind(se.matrix,se.vector)
			abcissa.matrix = cbind(abcissa.matrix,abcissa.vector)
		}
	}
	error.plus=mean.matrix+se.matrix
	error.minus=mean.matrix-se.matrix
	top=max(mean.matrix)
	top.se=max(se.matrix)
	ylim=c(0,top+3*top.se)
	names.arg=unique(as.factor(variable2))
	barplot(mean.matrix,beside=T,xlim=xlim,ylim=ylim,names.arg=names.arg,legend.text=legend.text,args.legend=args.legend,col=col,main=main,sub=sub,xlab=xlab,ylab=ylab)
	errbar(abcissa.matrix,mean.matrix,error.plus,error.minus,add=T)
}

##### SCATTER PLOT WITH TREND LINE
scatter.and.trend <- function (variable1, variable2, main=NULL, sub=NULL, xlab="", ylab="") {
	plot(variable1, variable2, main=main, sub=sub, xlab=xlab, ylab=ylab)
	model=lm(variable2~variable1)
	intercept=model$coefficients[1]
	slope=model$coefficients[2]
	abline(intercept,slope)
}