R/gdm.variable.importance.R
In fitzLab-AL/gdm: Generalized Dissimilarity Modeling
Defines functions
gdm.varImp
Documented in
gdm.varImp
#' @title Assess Predictor Importance and Quantify Model Significance in a Fitted
#' Generalized Dissimilarity Model.
#'
#' @description This function uses matrix permutation to perform model and
#' predictor significance testing and to estimate predictor importance in a
#' generalized dissimilarity model. The function can be run in parallel on
#' multicore machines to reduce computation time.
#'
#' @usage gdm.varImp(spTable, geo, splines = NULL, knots = NULL,
#' predSelect = FALSE, nPerm = 50, pValue=0.05, parallel = FALSE, cores = 2,
#' sampleSites = 1, sampleSitePairs = 1, outFile = NULL)
#'
#' @param spTable A site-pair table, same as used to fit a \code{\link[gdm]{gdm}}.
#'
#' @param geo Similar to the \code{\link[gdm]{gdm}} geo argument. The only
#' difference is that the geo argument does not have a default in this function.
#'
#' @param splines Same as the \code{\link[gdm]{gdm}} splines argument. Note that
#' the current implementation requires that all predictors have the same number of
#' splines.
#'
#' @param knots Same as the \code{\link[gdm]{gdm}} knots argument.
#'
#' @param predSelect Set to TRUE to perform predictor selection using matrix
#' permutation and backward elimination. Default is FALSE. When predSelect = FALSE
#' results will be returned only for a model fit with all predictors.
#'
#' @param nPerm Number of permutations to use to estimate p-values. Default is 50.
#'
#' @param pValue The p-value to use for predictor selection / elimination. Default is 0.05.
#'
#' @param parallel Whether or not to run the matrix permutations and model
#' fitting in parallel. Parallel processing is highly recommended when either
#' (i) the nPerms argument is large (>100) or (ii) a large number of site-pairs
#' (and / or variables) are being used in model fitting (note computation demand
#' can be reduced using subsampling - see next arguments). The default is FALSE.
#'
#' @param cores When the parallel argument is set to TRUE, the number of cores
#' to be registered for parallel processing. Must be <= the number of cores in
#' the machine running the function. There is no benefit to setting the number of
#' cores greater than the number of predictors in the model.
#'
#' @param sampleSites The fraction (0-1, though a value of 0 would be silly,
#' wouldn't it?) of \emph{sites to retain} from the full site-pair table. If
#' less than 1, this argument will completely remove a fraction of sites such
#' that they are not used in the permutation routines.
#'
#' @param sampleSitePairs The fraction (0-1) of \emph{site-pairs (i.e., rows)
#' to retain} from the full site-pair table - in other words, all sites will
#' be used in the permutation routines (assuming sampleSites = 1), but not
#' all \emph{site-pair combinations}. In the case where both the sampleSites
#' and the sampleSitePairs argument have values less than 1, sites first will
#' be removed using the sampleSites argument, followed by removal of site-pairs
#' using the sampleSitePairs argument. Note that the number of site-pairs
#' removed is based on the fraction of the resulting site-pair table after
#' sites have been removed, not on the size of the full site-pair table.
#'
#' @param outFile An optional character string to write the object returned by
#' the function to disk as an .RData object (".RData" is not required as part
#' of the file name). The .RData object will contain a single list with the
#' name of "outObject". The default is NULL, meaning that no file will be written.
#'
#' @details To test model significance, first a model is fit using all predictors and
#' un-permuted environmental data. Any predictor for which the sum of the I-spline
#' coefficients sum to zero is preemptively removed. Next, the environmental data are permuted
#' nPerm times (by randomizing the order of the rows) and a GDM is fit to each
#' permuted table. Model significance is determined by comparing the deviance
#' explained by GDM fit to the un-permuted table to the distribution of deviance
#' explained values from GDM fit to the nPerm permuted tables. To assess predictor
#' significance, this process is repeated for each predictor individually (i.e.,
#' only the data for the predictor being tested is permuted rather than the entire
#' environmental table). Predictor importance is quantified as the percent change
#' in deviance explained between a model fit with and without that predictor permuted. If
#' predSelect=TRUE, this process continues by next permutating the site-pair
#' table nPerm times, but removing one predictor at a time and reassessing
#' predictor importance and significance. At each step, the least important
#' predictor is dropped (backward elimination) and the process continues until
#' all non-significant predictors are removed, with significance level being set
#' by the user and the pValue argument.
#'
#' @return A list of four tables. The first table summarizes full model deviance,
#' percent deviance explained by the full model, the p-value of the full model,
#' and the number of permutations used to calculate the statistics for each
#' fitted model (i.e., the full model and each model with predictors removed in
#' succession during the backward elimination procedure if predSelect=T). The
#' remaining three tables summarize (1) predictor importance, (2) predictor
#' significance, and (3) the number of permutations used to calculate the
#' statistics for that model, which is provided because some GDMs may fail
#' to converge for some permutations / predictor combinations and you might want to
#' know how many permutations were used when calculating statistics. Or maybe you
#' don't, you decide.
#'
#' Predictor importance is measured as the percent decrease in deviance explained
#' between the full model and the deviance explained by a model fit with that predictor
#' permuted. Significance is estimated using the bootstrapped p-value when the
#' predictor has been permuted. For most cases, the number of permutations will
#' equal the nPerm argument. However, the value may be less should any of the models
#' fit to them permuted tables fail to converge.
#'
#' If predSelect=FALSE, the tables will have values only in the first column.
#'
#' @author Matt Fitzpatrick and Karel Mokany
#'
#' @references Ferrier S, Manion G, Elith J, Richardson, K (2007) Using
#' generalized dissimilarity modelling to analyse and predict patterns of
#' beta diversity in regional biodiversity assessment. \emph{Diversity &
#' Distributions} 13, 252-264.
#'
#' Fitzpatrick, MC, Sanders NJ, Ferrier S, Longino JT, Weiser MD, and RR Dunn.
#' 2011. Forecasting the Future of Biodiversity: a Test of Single- and
#' Multi-Species Models for Ants in North America. \emph{Ecography} 34: 836-47.
#'
#' @examples
#' ##fit table environmental data
#' ##set up site-pair table using the southwest data set
#' sppData <- southwest[, c(1,2,13,14)]
#' envTab <- southwest[, c(2:ncol(southwest))]
#' sitePairTab <- formatsitepair(sppData, 2, XColumn="Long", YColumn="Lat",
#' sppColumn="species", siteColumn="site", predData=envTab)
#'
#' ## not run
#' #modTest <- gdm.varImp(sitePairTab, geo=T, nPerm=50, parallel=T, cores=10, predSelect=T)
#' #barplot(sort(modTest$`Predictor Importance`[,1], decreasing=T))
#'
#' @keywords gdm
#'
#' @importFrom pbapply pbreplicate
#' @importFrom parallel makeCluster
#' @importFrom parallel stopCluster
#' @importFrom foreach %dopar%
#' @importFrom foreach foreach
#' @importFrom doParallel registerDoParallel
#' @importFrom utils txtProgressBar
#' @importFrom utils setTxtProgressBar
#'
#' @export
gdm.varImp <- function(spTable, geo, splines=NULL, knots=NULL, predSelect=FALSE,
nPerm=50, pValue=0.05, parallel=FALSE, cores=2, sampleSites=1,
sampleSitePairs=1, outFile=NULL){
##assign k to prevent issues with cran checking
k <- NULL
##error checking for input objects
##checks to see if in site-pair format from formatsitepair function
if(!is(spTable, "gdmData")){
warning("The spTable object is not of class 'gdmData'. See the formatsitepair function for help.")
}
##checks to makes sure data is a matrix or data frame
if(!(is(spTable, "gdmData") | is(spTable, "matrix") | is(spTable, "data.frame"))){
stop("spTable argument needs to be of class 'gdmData', 'matrix', or 'data frame'")
}
##sanity check on the data table
if(ncol(spTable) < 6){
stop("spTable object requires at least 6 columns: distance, weights, s1.xCoord, s1.yCoord, s2.xCoord, s2.yCoord")
}
if(nrow(spTable) < 1){
stop("The spTable object contains zero rows of data.")
}
##checks that geo has either TRUE or FALSE
if(!(geo==TRUE | geo==FALSE)){
stop("The geo argument must be either TRUE or FALSE.")
}
##makes sure splines is a numeric vector
if(is.null(splines)==FALSE & !is(splines, "numeric")){
stop("The splines argument needs to be a numeric data type.")
}
##checks knots inputs
if(is.null(knots)==FALSE & !is(knots, "numeric")){
stop("The knots argument needs to be a numeric data type.")
}
##checks that predSelect has either TRUE or FALSE
if(!(predSelect==TRUE | predSelect==FALSE)){
stop("The predSelect argument must be either TRUE or FALSE.")
}
##makes sure that nPerm is a positive integer
if((is.null(nPerm)==FALSE & is.numeric(nPerm)==FALSE) | nPerm<1){
stop("The nPerm argument needs to be a positive integer.")
}
##checks that parallel has either TRUE or FALSE
if(!(parallel==TRUE | parallel==FALSE)){
stop("The parallel argument must be either TRUE or FALSE.")
}
##makes sure that cores has a value when parallel is true
if(parallel==TRUE & is.null(cores)==TRUE){
stop("If parallel==TRUE, the number of cores must be specified.")
}
##makes sure that cores is a positive integer
if((is.null(cores)==FALSE & is.numeric(cores)==FALSE) | cores<1){
stop("The cores argument needs to be a positive integer.")
}
##makes sure that both sampleSites and sampleSitePairs are a number between 0 and 1,
##and that neither is equal to 0
if(is.numeric(sampleSites)==FALSE | sampleSites<0 | sampleSites>1){
stop("The sampleSites argument needs to be a positive number between 0 and 1.")
}
if(is.numeric(sampleSitePairs)==FALSE | sampleSitePairs<0 | sampleSitePairs>1){
stop("The sampleSitePairs argument needs to be a positive number between 0 and 1.")
}
if(sampleSites==0){
stop("A sampleSites value of 0 will remove all sites from the analysis.")
}
if(sampleSitePairs==0){
stop("A sampleSitePairs value of 0 will remove all sites from the analysis.")
}
##checks to see if the user has requested for an output file to be written, and if so
##makes sure that it is formatted correctly
if(is.null(outFile)==FALSE){
##first makes sure outFile is a string
if(is.character(outFile)==FALSE){
stop("The outFile argument needs to be a character string of the directory and file name you wish the tables to be written to")
}
##makes sure that text has ".RData" in it, if not, adds it
outFileChar <- nchar(outFile)
if(substr(outFile, outFileChar-5, outFileChar)!=".RData"){
outFile <- paste(outFile, ".RData", sep="")
}
##checks to see if there is a path as well as a file name
if(length(strsplit(outFile,"/")[[1]])>1){
splitOutFile <- strsplit(outFile,"/")[[1]][-length(strsplit(outFile,"/")[[1]])]
dir.create(paste(splitOutFile, collapse="/"))
}else{
outFile <- paste("./", outFile, sep="")
}
}
##double makes sure these values are integers, seems to truncate if not
nPerm <- as.integer(nPerm)
cores <- as.integer(cores)
##removes a user specified number of sites from the site-pair table
if(sampleSites<1){
spTable <- subsample.sitepair(spTable, sampleSites=sampleSites)
##throws warning if sampleSitePairs<1 as well
if(sampleSitePairs<1){
warning("You have selected to randomly remove sites and/or site-pairs.")
}
}
##removes a user specified number of site-pairs from the site-pair table
if(sampleSitePairs<1){
##determine which rows to remove
numRm <- sample(1:nrow(spTable), round(nrow(spTable)*(1-sampleSitePairs)))
spTable <- spTable[-c(numRm),]
}
##check that the response data is [0..1]
rtmp <- spTable[,1]
if(length(rtmp[rtmp<0]) > 0){
stop("The spTable contains negative distance values. Must be between 0 - 1.")
}
if (length(rtmp[rtmp>1]) > 0){
stop("The spTable contains distance values greater than 1. Must be between 0 - 1.")
}
# number of variables in the site-pair table, adds 1 if geo=TRUE
nVars <- (ncol(spTable)-6)/2
# create vector of variable names
varNames <- colnames(spTable[c(7:(6+nVars))])
varNames <- sapply(strsplit(varNames, "s1\\."), "[[", 2)
if(geo==TRUE){
nVars <- nVars + 1
varNames <- c("Geographic", varNames)
}
if(nVars<2){
stop("Function requires at least two predictor variables.")
}
# run initial GDM to see if any vars have zero I-spline coeffs
message(paste0("Fitting initial model with all ", nVars, " predictors..."))
Sys.sleep(0.5)
fullGDM <- gdm(spTable, geo=geo, splines=splines, knots=knots)
# check for zero coeffs
thiscoeff <- 1
thisquant <- 1
sumCoeff <- NULL
for(i in 1:length(fullGDM$predictors)){
numsplines <- fullGDM$splines[[i]]
holdCoeff <- NULL
for(j in 1:numsplines){
holdCoeff[j] <- fullGDM$coefficients[[thiscoeff]]
thiscoeff <- thiscoeff + 1
}
sumCoeff[i] <- sum(holdCoeff)
}
# remove any predictors with sumCoeff=0
zeroSum <- fullGDM$predictors[which(sumCoeff==0)]
if(length(zeroSum)>0){
for(p in 1:length(zeroSum)){
message(paste0("Sum of I-spline coefficients for predictor ", zeroSum[p]," = 0"))
Sys.sleep(0.5)}
#message("\n")
for(p in 1:length(zeroSum)){
if(zeroSum[p]=="Geographic"){
message("Setting Geo=FALSE and proceeding with permutation testing...")
} else {
message(paste0("Removing ", zeroSum[p], " and proceeding with permutation testing..."))
}
Sys.sleep(0.5)
}
if(length(grep("Geographic", zeroSum))==1){
geo <- FALSE
zeroSum <- zeroSum[-grep("Geographic", zeroSum)]
}
for(z in zeroSum){
##select variable columns to be removed from original site-pair table
testVarCols1 <- grep(paste("^s1\\.", z, "$", sep=""), colnames(spTable))
testVarCols2 <- grep(paste("^s2\\.", z, "$", sep=""), colnames(spTable))
spTable <- spTable[,-c(testVarCols1, testVarCols2)]
}
}
# recalculate the number of variables in the site-pair table, adds 1 if geo=TRUE
nVars <- (ncol(spTable)-6)/2
# create vector of variable names
varNames <- colnames(spTable[c(7:(6+nVars))])
varNames <- sapply(strsplit(varNames, "s1\\."), "[[", 2)
if(geo==TRUE){
nVars <- nVars + 1
varNames <- c("Geographic", varNames)
}
# stop routine if fewer than two variables remain
if(nVars<2){
stop("Function requires at least two predictor variables.")
}
# reduce number of cores to nVars
if(cores>nVars){
cores <- nVars
}
# crate spline object
splines <- rep(unique(splines), nVars)
##First create a spTable to determine the index of each site in the site-pair table
sortMatX <- sapply(1:nrow(spTable), function(i, spTab){c(spTab[i,3], spTab[i,5])}, spTab=spTable)
sortMatY <- sapply(1:nrow(spTable), function(i, spTab){c(spTab[i,4], spTab[i,6])}, spTab=spTable)
sortMatNum <- sapply(1:nrow(spTable), function(i){c(1,2)})
sortMatRow <- sapply(1:nrow(spTable), function(i){c(i,i)})
##adds a column of NA for index to be added to
fullSortMat <- cbind(as.vector(sortMatX), as.vector(sortMatY), as.vector(sortMatNum), as.vector(sortMatRow), rep(NA, length(sortMatX)))
##assigns sites by unique coordinates
siteByCoords <- as.data.frame(unique(fullSortMat[,1:2]))
##number of sites to expect by coordinates
numSites <- nrow(siteByCoords)
##assigns site index based on coordinates
for(i in 1:numSites){
fullSortMat[which(fullSortMat[,1]==siteByCoords[i,1] & fullSortMat[,2]==siteByCoords[i,2]),5] <- i
}
##create index table to know where each site is in input site-pair table
indexTab <- matrix(NA,nrow(spTable),2)
for(iRow in 1:nrow(fullSortMat)){
indexTab[fullSortMat[iRow,4],fullSortMat[iRow,3]] <- fullSortMat[iRow,5]
}
## And remove the sorting table and supporting objects to free up memory
rm(fullSortMat)
rm(sortMatX)
rm(sortMatY)
rm(sortMatNum)
rm(sortMatRow)
rm(siteByCoords)
##create site x predictor table, to be able to rebuild site-pair table later in function
exBySite <- lapply(1:numSites, function(i, index, tab){
rowSites <- which(index[,1] %in% i)
if(length(rowSites)<1){
rowSites <- which(index[,2] %in% i)
}
exSiteData <- tab[rowSites[1],]
return(exSiteData)
}, index=indexTab, tab=spTable)
##identifies the one site not in the first column of the index table
outSite <- which(!(1:numSites %in% indexTab[,1]))
##sets up siteXvar table, uses for loop to make sure have steps correct
for(i in 1:length(exBySite)){
##grabs row and identify if should take s1 or s2 by rather or not number appears in outsite
siteRow <- exBySite[[i]]
if(i %in% outSite){
##extracts the data from the site-pair table by site
siteRow <- siteRow[grep("s2\\.", colnames(siteRow))]
colnames(siteRow) <- sapply(strsplit(colnames(siteRow), "s2\\."), "[[", 2)
}else{
##extracts the data from the site-pair table by site
siteRow <- siteRow[grep("s1\\.", colnames(siteRow))]
colnames(siteRow) <- sapply(strsplit(colnames(siteRow), "s1\\."), "[[", 2)
}
exBySite[[i]] <- siteRow
}
##transforms data from list to table
siteData <- do.call("rbind", exBySite)
##sets up objects to be returned by the function
modelTestValues <- matrix(NA,4,nVars*5,dimnames = list(c("Model deviance",
"Percent deviance explained",
"Model p-value",
"Fitted permutations"),
c("All predictors",
"1-removed",
paste(seq(2,nVars*5-1), "-removed", sep=""))))
##deviance reduction predictor table
varImpTable <- matrix(NA, nVars, nVars*5-1)
rownames(varImpTable) <- varNames
colnames(varImpTable) <- c("All predictors",
"1-removed",
paste(seq(2,nVars*5-2), "-removed", sep=""))
##p value predictor table
pValues <- nModsConverge <- varImpTable
##assigns given site-pair table to new variable, to prevent changing the original input
currSitePair <- spTable
nullGDMFullFit <- 0 ##a variable to track if the fully fitted gdm model returned a NULL object
# create the set up permuted site-pair tables to be used for all
# downstream analyses
message(paste0("Creating ", nPerm, " permuted site-pair tables..."))
# use replicate if number of perms is small or parallel == FALSE
if(parallel==F | nPerm <= 25){
permSpt <- pbreplicate(nPerm, list(permutateSitePair(currSitePair,
siteData,
indexTab,
varNames)))}
# create permuted tables in parallel otherwise
if(parallel==T & nPerm > 25){
# set up parallel processing
cl <- makeCluster(cores)
registerDoParallel(cl)
permSpt <- foreach(k=1:nPerm,
.verbose=F,
.packages=c("gdm"),
.export=c("permutateSitePair")) %dopar%
permutateSitePair(currSitePair, siteData, indexTab, varNames)
stopCluster(cl)}
# create new vector to track varNames
varNames.x <- varNames
message("Starting model assessment...")
for(v in 1:length(varNames)){
# ends the loop if only 1 variable remains
if(length(varNames.x)<2){
stop("Only one predictor remains...variable assessment stopped.")
}
if(v>1){
message(paste0("Removing ", names(elimVar), " and proceeding with the next round of permutations."))
}
if(is.numeric(splines)){
splines <- rep(unique(splines), length(varNames.x))
}
# runs gdm, first time on site-pair table with all variables
# As variables are eliminated the "full" site-pair table will
# contain fewer variables
fullGDM <- gdm(currSitePair, geo=geo, splines=splines, knots=knots)
message(paste0("Percent deviance explained by the full model = ", round(fullGDM$explained,3)))
if(is.null(fullGDM)==TRUE){
warning(paste("The model did not converge when testing variable: ", varNames.x[v],
". Terminating analysis and returning output completed up to this point.", sep=""))
break
}
# fit gdm to each permuted site-pair table
message("Fitting GDMs to the permuted site-pair tables...")
permGDM <- lapply(permSpt, function(x){
gdm(x, geo=geo, splines=splines, knots=knots)
})
##extracts deviance of permuted gdms
permModelDev <- sapply(permGDM, function(mod){mod$gdmdeviance})
##if needed, removes nulls from output permModelDev
modPerms <- length(which(sapply(permModelDev,is.null)==TRUE))
if(modPerms>0){
permModelDev <- unlist(permModelDev[-(which(sapply(permModelDev,is.null)==T))])
}
if(v>1){
colnames(modelTestValues)[v] <- colnames(pValues)[v] <- colnames(varImpTable)[v] <- colnames(nModsConverge)[v] <- paste0(names(elimVar),"-removed")
}
##begins to fill in the output table with data from fully fitted model
modelTestValues[1,v] <- round(fullGDM$gdmdeviance,3)
modelTestValues[2,v] <- round(fullGDM$explained,3)
#p-value
modelTestValues[3,v] <- round(sum(permModelDev<=fullGDM$gdmdeviance)/(nPerm-modPerms),3)
##fitted permutations
modelTestValues[4,v] <- round(nPerm-modPerms,0)
# now permutate each variable individually and fit models
if(parallel==TRUE){
if(length(varNames.x)<cores){
cores <- length(varNames.x)
}
# set up parallel processing
cl <- makeCluster(cores)
registerDoParallel(cl)
# foreach function to create site-pair tables with each variable permuted,
# fit gdms and extract deviance.
permVarDev <- foreach(k=1:length(varNames.x),
.verbose=F,
.packages=c("gdm"),
.export = c("currSitePair")) %dopar%{
if(varNames.x[k]!="Geographic"){
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
idx1 <- grep(paste("^s1\\.", varNames.x[k], "$", sep=""), colnames(x))
idx2 <- grep(paste("^s2\\.", varNames.x[k], "$", sep=""), colnames(x))
spt[,c(idx1, idx2)] <- x[,c(idx1, idx2)]
return(spt)
})}
if(varNames.x[k]=="Geographic"){
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
s1 <- sample(1:nrow(spt), nrow(spt))
s2 <- sample(1:nrow(spt), nrow(spt))
s3 <- sample(1:nrow(spt), nrow(spt))
s4 <- sample(1:nrow(spt), nrow(spt))
spt[,3] <- spt[s1,3]
spt[,4] <- spt[s2,4]
spt[,5] <- spt[s3,5]
spt[,6] <- spt[s4,6]
return(spt)
})}
gdmPermVar <- lapply(lll, function(x){
try(gdm(x, geo=geo, splines=splines, knots=knots))
})
##extracts deviance of permuted gdms
permModelDev <- sapply(gdmPermVar, function(mod){mod$gdmdeviance})
return(permModelDev)
}
##closes cores
#close(pb)
stopCluster(cl)
}
if(parallel==FALSE){
# for-loop to create site-pair tables with each variable permuted,
# fit gdms and extract deviance.
permVarDev <- list()
for(k in 1:length(varNames.x)){
if(varNames.x[k]!="Geographic"){
message(paste0("Assessing importance of ", varNames.x[k], "..."))
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
idx1 <- grep(paste("^s1\\.", varNames.x[k], "$", sep=""), colnames(x))
idx2 <- grep(paste("^s2\\.", varNames.x[k], "$", sep=""), colnames(x))
spt[,c(idx1, idx2)] <- x[,c(idx1, idx2)]
return(spt)
})}
if(varNames.x[k]=="Geographic"){
message("Assessing importance of geographic distance...")
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
s1 <- sample(1:nrow(spt), nrow(spt))
s2 <- sample(1:nrow(spt), nrow(spt))
s3 <- sample(1:nrow(spt), nrow(spt))
s4 <- sample(1:nrow(spt), nrow(spt))
spt[,3] <- spt[s1,3]
spt[,4] <- spt[s2,4]
spt[,5] <- spt[s3,5]
spt[,6] <- spt[s4,6]
return(spt)
})}
gdmPermVar <- lapply(lll, function(x){
try(gdm(x, geo=geo, splines=splines, knots=knots))
})
##extracts deviance of permuted gdms
permVarDev[[k]] <- unlist(sapply(gdmPermVar, function(mod){mod$gdmdeviance}))
}
}
names(permVarDev) <- varNames.x
nullDev <- fullGDM$nulldeviance
for(var in varNames.x){
grepper <- grep(paste0("^",var,"$"), names(permVarDev))
varDevTab <- permVarDev[[grepper]]
# number of perms for which GDM converged
nConv <- length(varDevTab)
nModsConverge[which(rownames(varImpTable) == var),v] <- nConv
# remove NULLs (GDMs that did not converge)
#if(nConv>0){
#varDevTab <- unlist(varDevTab[-(which(sapply(is.null(varDevTab))))])
#}
# calculate variable importance
varDevExplained <- 100*(nullDev-varDevTab)/nullDev
varImpTable[which(rownames(varImpTable) == var),v] <- median(100 * abs((varDevExplained - fullGDM$explained)/fullGDM$explained))
if(var!="Geographic"){
##select variable columns to be removed from original site-pair table
testVarCols1 <- grep(paste("^s1\\.", var, "$", sep=""), colnames(currSitePair))
testVarCols2 <- grep(paste("^s2\\.", var, "$", sep=""), colnames(currSitePair))
testSitePair <- currSitePair[,-c(testVarCols1, testVarCols2)]
##run gdm for the missing variable
noVarGDM <- gdm(testSitePair, geo=geo, splines=splines[-1], knots=knots)
} else {
noVarGDM <- gdm(currSitePair, geo=F, splines=splines[-1], knots=knots)
}
# calculate p-Value
permDevReduct <- noVarGDM$gdmdeviance - varDevTab
pValues[which(rownames(pValues) == var),v] <- sum(permDevReduct>=(varDevTab - fullGDM$gdmdeviance))/(nConv)
}
if(max(na.omit(pValues[,v]))<pValue){
message("All remaining predictors are significant, ceasing assessment.")
message(paste0("Percent deviance explained by final model = ", round(fullGDM$explained,3)))
message("Final set of predictors returned: ")
for(vvv in 1:length(fullGDM$predictors)){
message(fullGDM$predictors[vvv])
}
break
}
if(predSelect==T){
# eliminate least important variable
#elimVar <- which.min(varImpTable[,v])
# eliminate variable with highest p-value
elimVar <- which.max(pValues[,v])
if(names(elimVar)!="Geographic"){
##select variable columns to be removed from original site-pair table
remVarCols1 <- grep(paste("^s1\\.", names(elimVar), "$", sep=""), colnames(currSitePair))
remVarCols2 <- grep(paste("^s2\\.", names(elimVar), "$", sep=""), colnames(currSitePair))
# remove columns from site-pair table
currSitePair <- currSitePair[,-c(remVarCols1, remVarCols2)]
# remove columns from permuted site-pair tables
permSpt <- lapply(permSpt, function(x){
x[,-c(remVarCols1, remVarCols2)]
})
} else {
geo <- F
}
varNames.x <- varNames.x[-which(varNames.x==names(elimVar))]
}
if(v==1 & predSelect==F){
message("Backwards elimination not selected by user (predSelect=F). Ceasing assessment.")
message(paste0("Percent deviance explained by final model = ", round(fullGDM$explained,3)))
message("Final set of predictors returned: ")
for(vvv in 1:length(fullGDM$predictors)){
message(fullGDM$predictors[vvv])
}
break
}
}
if(v==1 & predSelect==F){
# Model assessment
modelTestVals <- data.frame(matrix(round(modelTestValues[,1], 3), ncol=1))
rownames(modelTestVals) <- rownames(modelTestValues)
colnames(modelTestVals) <- "All predictors"
#Variable importance
varImpTab <- data.frame(matrix(round(varImpTable[,1], 3), ncol=1))
rownames(varImpTab) <- rownames(varImpTable)
colnames(varImpTab) <- "All predictors"
#Variable selection p-values
pVals <- varImpTab
pVals[,1] <- round(pValues[,1], 3)
#Variable selection model convergence
nModsConv <- varImpTab
nModsConv[,1] <- round(nModsConverge[,1], 3)
outObject <- list(modelTestVals, varImpTab, pVals, nModsConv)
names(outObject) <- c("Model assessment", "Predictor Importance", "Predictor p-values", "Model Convergence")
} else{
outObject <- list(round(modelTestValues[,1:v], 3), round(varImpTable[,1:v],3), round(pValues[,1:v],3), nModsConverge[,1:v])
names(outObject) <- c("Model assessment", "Predictor Importance", "Predictor p-values", "Model Convergence")
}
##if given, writes out files to space on disk
if(is.null(outFile)==FALSE){
save(outObject, file=outFile)
}
return(outObject)
}
fitzLab-AL/gdm documentation built on April 17, 2025, 11 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions gdm.varImp
Documented in gdm.varImp
#' @title Assess Predictor Importance and Quantify Model Significance in a Fitted
#' Generalized Dissimilarity Model.
#'
#' @description This function uses matrix permutation to perform model and
#' predictor significance testing and to estimate predictor importance in a
#' generalized dissimilarity model. The function can be run in parallel on
#' multicore machines to reduce computation time.
#'
#' @usage gdm.varImp(spTable, geo, splines = NULL, knots = NULL,
#' predSelect = FALSE, nPerm = 50, pValue=0.05, parallel = FALSE, cores = 2,
#' sampleSites = 1, sampleSitePairs = 1, outFile = NULL)
#'
#' @param spTable A site-pair table, same as used to fit a \code{\link[gdm]{gdm}}.
#'
#' @param geo Similar to the \code{\link[gdm]{gdm}} geo argument. The only
#' difference is that the geo argument does not have a default in this function.
#'
#' @param splines Same as the \code{\link[gdm]{gdm}} splines argument. Note that
#' the current implementation requires that all predictors have the same number of
#' splines.
#'
#' @param knots Same as the \code{\link[gdm]{gdm}} knots argument.
#'
#' @param predSelect Set to TRUE to perform predictor selection using matrix
#' permutation and backward elimination. Default is FALSE. When predSelect = FALSE
#' results will be returned only for a model fit with all predictors.
#'
#' @param nPerm Number of permutations to use to estimate p-values. Default is 50.
#'
#' @param pValue The p-value to use for predictor selection / elimination. Default is 0.05.
#'
#' @param parallel Whether or not to run the matrix permutations and model
#' fitting in parallel. Parallel processing is highly recommended when either
#' (i) the nPerms argument is large (>100) or (ii) a large number of site-pairs
#' (and / or variables) are being used in model fitting (note computation demand
#' can be reduced using subsampling - see next arguments). The default is FALSE.
#'
#' @param cores When the parallel argument is set to TRUE, the number of cores
#' to be registered for parallel processing. Must be <= the number of cores in
#' the machine running the function. There is no benefit to setting the number of
#' cores greater than the number of predictors in the model.
#'
#' @param sampleSites The fraction (0-1, though a value of 0 would be silly,
#' wouldn't it?) of \emph{sites to retain} from the full site-pair table. If
#' less than 1, this argument will completely remove a fraction of sites such
#' that they are not used in the permutation routines.
#'
#' @param sampleSitePairs The fraction (0-1) of \emph{site-pairs (i.e., rows)
#' to retain} from the full site-pair table - in other words, all sites will
#' be used in the permutation routines (assuming sampleSites = 1), but not
#' all \emph{site-pair combinations}. In the case where both the sampleSites
#' and the sampleSitePairs argument have values less than 1, sites first will
#' be removed using the sampleSites argument, followed by removal of site-pairs
#' using the sampleSitePairs argument. Note that the number of site-pairs
#' removed is based on the fraction of the resulting site-pair table after
#' sites have been removed, not on the size of the full site-pair table.
#'
#' @param outFile An optional character string to write the object returned by
#' the function to disk as an .RData object (".RData" is not required as part
#' of the file name). The .RData object will contain a single list with the
#' name of "outObject". The default is NULL, meaning that no file will be written.
#'
#' @details To test model significance, first a model is fit using all predictors and
#' un-permuted environmental data. Any predictor for which the sum of the I-spline
#' coefficients sum to zero is preemptively removed. Next, the environmental data are permuted
#' nPerm times (by randomizing the order of the rows) and a GDM is fit to each
#' permuted table. Model significance is determined by comparing the deviance
#' explained by GDM fit to the un-permuted table to the distribution of deviance
#' explained values from GDM fit to the nPerm permuted tables. To assess predictor
#' significance, this process is repeated for each predictor individually (i.e.,
#' only the data for the predictor being tested is permuted rather than the entire
#' environmental table). Predictor importance is quantified as the percent change
#' in deviance explained between a model fit with and without that predictor permuted. If
#' predSelect=TRUE, this process continues by next permutating the site-pair
#' table nPerm times, but removing one predictor at a time and reassessing
#' predictor importance and significance. At each step, the least important
#' predictor is dropped (backward elimination) and the process continues until
#' all non-significant predictors are removed, with significance level being set
#' by the user and the pValue argument.
#'
#' @return A list of four tables. The first table summarizes full model deviance,
#' percent deviance explained by the full model, the p-value of the full model,
#' and the number of permutations used to calculate the statistics for each
#' fitted model (i.e., the full model and each model with predictors removed in
#' succession during the backward elimination procedure if predSelect=T). The
#' remaining three tables summarize (1) predictor importance, (2) predictor
#' significance, and (3) the number of permutations used to calculate the
#' statistics for that model, which is provided because some GDMs may fail
#' to converge for some permutations / predictor combinations and you might want to
#' know how many permutations were used when calculating statistics. Or maybe you
#' don't, you decide.
#'
#' Predictor importance is measured as the percent decrease in deviance explained
#' between the full model and the deviance explained by a model fit with that predictor
#' permuted. Significance is estimated using the bootstrapped p-value when the
#' predictor has been permuted. For most cases, the number of permutations will
#' equal the nPerm argument. However, the value may be less should any of the models
#' fit to them permuted tables fail to converge.
#'
#' If predSelect=FALSE, the tables will have values only in the first column.
#'
#' @author Matt Fitzpatrick and Karel Mokany
#'
#' @references Ferrier S, Manion G, Elith J, Richardson, K (2007) Using
#' generalized dissimilarity modelling to analyse and predict patterns of
#' beta diversity in regional biodiversity assessment. \emph{Diversity &
#' Distributions} 13, 252-264.
#'
#' Fitzpatrick, MC, Sanders NJ, Ferrier S, Longino JT, Weiser MD, and RR Dunn.
#' 2011. Forecasting the Future of Biodiversity: a Test of Single- and
#' Multi-Species Models for Ants in North America. \emph{Ecography} 34: 836-47.
#'
#' @examples
#' ##fit table environmental data
#' ##set up site-pair table using the southwest data set
#' sppData <- southwest[, c(1,2,13,14)]
#' envTab <- southwest[, c(2:ncol(southwest))]
#' sitePairTab <- formatsitepair(sppData, 2, XColumn="Long", YColumn="Lat",
#' sppColumn="species", siteColumn="site", predData=envTab)
#'
#' ## not run
#' #modTest <- gdm.varImp(sitePairTab, geo=T, nPerm=50, parallel=T, cores=10, predSelect=T)
#' #barplot(sort(modTest$`Predictor Importance`[,1], decreasing=T))
#'
#' @keywords gdm
#'
#' @importFrom pbapply pbreplicate
#' @importFrom parallel makeCluster
#' @importFrom parallel stopCluster
#' @importFrom foreach %dopar%
#' @importFrom foreach foreach
#' @importFrom doParallel registerDoParallel
#' @importFrom utils txtProgressBar
#' @importFrom utils setTxtProgressBar
#'
#' @export
gdm.varImp <- function(spTable, geo, splines=NULL, knots=NULL, predSelect=FALSE,
nPerm=50, pValue=0.05, parallel=FALSE, cores=2, sampleSites=1,
sampleSitePairs=1, outFile=NULL){
##assign k to prevent issues with cran checking
k <- NULL
##error checking for input objects
##checks to see if in site-pair format from formatsitepair function
if(!is(spTable, "gdmData")){
warning("The spTable object is not of class 'gdmData'. See the formatsitepair function for help.")
}
##checks to makes sure data is a matrix or data frame
if(!(is(spTable, "gdmData") | is(spTable, "matrix") | is(spTable, "data.frame"))){
stop("spTable argument needs to be of class 'gdmData', 'matrix', or 'data frame'")
}
##sanity check on the data table
if(ncol(spTable) < 6){
stop("spTable object requires at least 6 columns: distance, weights, s1.xCoord, s1.yCoord, s2.xCoord, s2.yCoord")
}
if(nrow(spTable) < 1){
stop("The spTable object contains zero rows of data.")
}
##checks that geo has either TRUE or FALSE
if(!(geo==TRUE | geo==FALSE)){
stop("The geo argument must be either TRUE or FALSE.")
}
##makes sure splines is a numeric vector
if(is.null(splines)==FALSE & !is(splines, "numeric")){
stop("The splines argument needs to be a numeric data type.")
}
##checks knots inputs
if(is.null(knots)==FALSE & !is(knots, "numeric")){
stop("The knots argument needs to be a numeric data type.")
}
##checks that predSelect has either TRUE or FALSE
if(!(predSelect==TRUE | predSelect==FALSE)){
stop("The predSelect argument must be either TRUE or FALSE.")
}
##makes sure that nPerm is a positive integer
if((is.null(nPerm)==FALSE & is.numeric(nPerm)==FALSE) | nPerm<1){
stop("The nPerm argument needs to be a positive integer.")
}
##checks that parallel has either TRUE or FALSE
if(!(parallel==TRUE | parallel==FALSE)){
stop("The parallel argument must be either TRUE or FALSE.")
}
##makes sure that cores has a value when parallel is true
if(parallel==TRUE & is.null(cores)==TRUE){
stop("If parallel==TRUE, the number of cores must be specified.")
}
##makes sure that cores is a positive integer
if((is.null(cores)==FALSE & is.numeric(cores)==FALSE) | cores<1){
stop("The cores argument needs to be a positive integer.")
}
##makes sure that both sampleSites and sampleSitePairs are a number between 0 and 1,
##and that neither is equal to 0
if(is.numeric(sampleSites)==FALSE | sampleSites<0 | sampleSites>1){
stop("The sampleSites argument needs to be a positive number between 0 and 1.")
}
if(is.numeric(sampleSitePairs)==FALSE | sampleSitePairs<0 | sampleSitePairs>1){
stop("The sampleSitePairs argument needs to be a positive number between 0 and 1.")
}
if(sampleSites==0){
stop("A sampleSites value of 0 will remove all sites from the analysis.")
}
if(sampleSitePairs==0){
stop("A sampleSitePairs value of 0 will remove all sites from the analysis.")
}
##checks to see if the user has requested for an output file to be written, and if so
##makes sure that it is formatted correctly
if(is.null(outFile)==FALSE){
##first makes sure outFile is a string
if(is.character(outFile)==FALSE){
stop("The outFile argument needs to be a character string of the directory and file name you wish the tables to be written to")
}
##makes sure that text has ".RData" in it, if not, adds it
outFileChar <- nchar(outFile)
if(substr(outFile, outFileChar-5, outFileChar)!=".RData"){
outFile <- paste(outFile, ".RData", sep="")
}
##checks to see if there is a path as well as a file name
if(length(strsplit(outFile,"/")[[1]])>1){
splitOutFile <- strsplit(outFile,"/")[[1]][-length(strsplit(outFile,"/")[[1]])]
dir.create(paste(splitOutFile, collapse="/"))
}else{
outFile <- paste("./", outFile, sep="")
}
}
##double makes sure these values are integers, seems to truncate if not
nPerm <- as.integer(nPerm)
cores <- as.integer(cores)
##removes a user specified number of sites from the site-pair table
if(sampleSites<1){
spTable <- subsample.sitepair(spTable, sampleSites=sampleSites)
##throws warning if sampleSitePairs<1 as well
if(sampleSitePairs<1){
warning("You have selected to randomly remove sites and/or site-pairs.")
}
}
##removes a user specified number of site-pairs from the site-pair table
if(sampleSitePairs<1){
##determine which rows to remove
numRm <- sample(1:nrow(spTable), round(nrow(spTable)*(1-sampleSitePairs)))
spTable <- spTable[-c(numRm),]
}
##check that the response data is [0..1]
rtmp <- spTable[,1]
if(length(rtmp[rtmp<0]) > 0){
stop("The spTable contains negative distance values. Must be between 0 - 1.")
}
if (length(rtmp[rtmp>1]) > 0){
stop("The spTable contains distance values greater than 1. Must be between 0 - 1.")
}
# number of variables in the site-pair table, adds 1 if geo=TRUE
nVars <- (ncol(spTable)-6)/2
# create vector of variable names
varNames <- colnames(spTable[c(7:(6+nVars))])
varNames <- sapply(strsplit(varNames, "s1\\."), "[[", 2)
if(geo==TRUE){
nVars <- nVars + 1
varNames <- c("Geographic", varNames)
}
if(nVars<2){
stop("Function requires at least two predictor variables.")
}
# run initial GDM to see if any vars have zero I-spline coeffs
message(paste0("Fitting initial model with all ", nVars, " predictors..."))
Sys.sleep(0.5)
fullGDM <- gdm(spTable, geo=geo, splines=splines, knots=knots)
# check for zero coeffs
thiscoeff <- 1
thisquant <- 1
sumCoeff <- NULL
for(i in 1:length(fullGDM$predictors)){
numsplines <- fullGDM$splines[[i]]
holdCoeff <- NULL
for(j in 1:numsplines){
holdCoeff[j] <- fullGDM$coefficients[[thiscoeff]]
thiscoeff <- thiscoeff + 1
}
sumCoeff[i] <- sum(holdCoeff)
}
# remove any predictors with sumCoeff=0
zeroSum <- fullGDM$predictors[which(sumCoeff==0)]
if(length(zeroSum)>0){
for(p in 1:length(zeroSum)){
message(paste0("Sum of I-spline coefficients for predictor ", zeroSum[p]," = 0"))
Sys.sleep(0.5)}
#message("\n")
for(p in 1:length(zeroSum)){
if(zeroSum[p]=="Geographic"){
message("Setting Geo=FALSE and proceeding with permutation testing...")
} else {
message(paste0("Removing ", zeroSum[p], " and proceeding with permutation testing..."))
}
Sys.sleep(0.5)
}
if(length(grep("Geographic", zeroSum))==1){
geo <- FALSE
zeroSum <- zeroSum[-grep("Geographic", zeroSum)]
}
for(z in zeroSum){
##select variable columns to be removed from original site-pair table
testVarCols1 <- grep(paste("^s1\\.", z, "$", sep=""), colnames(spTable))
testVarCols2 <- grep(paste("^s2\\.", z, "$", sep=""), colnames(spTable))
spTable <- spTable[,-c(testVarCols1, testVarCols2)]
}
}
# recalculate the number of variables in the site-pair table, adds 1 if geo=TRUE
nVars <- (ncol(spTable)-6)/2
# create vector of variable names
varNames <- colnames(spTable[c(7:(6+nVars))])
varNames <- sapply(strsplit(varNames, "s1\\."), "[[", 2)
if(geo==TRUE){
nVars <- nVars + 1
varNames <- c("Geographic", varNames)
}
# stop routine if fewer than two variables remain
if(nVars<2){
stop("Function requires at least two predictor variables.")
}
# reduce number of cores to nVars
if(cores>nVars){
cores <- nVars
}
# crate spline object
splines <- rep(unique(splines), nVars)
##First create a spTable to determine the index of each site in the site-pair table
sortMatX <- sapply(1:nrow(spTable), function(i, spTab){c(spTab[i,3], spTab[i,5])}, spTab=spTable)
sortMatY <- sapply(1:nrow(spTable), function(i, spTab){c(spTab[i,4], spTab[i,6])}, spTab=spTable)
sortMatNum <- sapply(1:nrow(spTable), function(i){c(1,2)})
sortMatRow <- sapply(1:nrow(spTable), function(i){c(i,i)})
##adds a column of NA for index to be added to
fullSortMat <- cbind(as.vector(sortMatX), as.vector(sortMatY), as.vector(sortMatNum), as.vector(sortMatRow), rep(NA, length(sortMatX)))
##assigns sites by unique coordinates
siteByCoords <- as.data.frame(unique(fullSortMat[,1:2]))
##number of sites to expect by coordinates
numSites <- nrow(siteByCoords)
##assigns site index based on coordinates
for(i in 1:numSites){
fullSortMat[which(fullSortMat[,1]==siteByCoords[i,1] & fullSortMat[,2]==siteByCoords[i,2]),5] <- i
}
##create index table to know where each site is in input site-pair table
indexTab <- matrix(NA,nrow(spTable),2)
for(iRow in 1:nrow(fullSortMat)){
indexTab[fullSortMat[iRow,4],fullSortMat[iRow,3]] <- fullSortMat[iRow,5]
}
## And remove the sorting table and supporting objects to free up memory
rm(fullSortMat)
rm(sortMatX)
rm(sortMatY)
rm(sortMatNum)
rm(sortMatRow)
rm(siteByCoords)
##create site x predictor table, to be able to rebuild site-pair table later in function
exBySite <- lapply(1:numSites, function(i, index, tab){
rowSites <- which(index[,1] %in% i)
if(length(rowSites)<1){
rowSites <- which(index[,2] %in% i)
}
exSiteData <- tab[rowSites[1],]
return(exSiteData)
}, index=indexTab, tab=spTable)
##identifies the one site not in the first column of the index table
outSite <- which(!(1:numSites %in% indexTab[,1]))
##sets up siteXvar table, uses for loop to make sure have steps correct
for(i in 1:length(exBySite)){
##grabs row and identify if should take s1 or s2 by rather or not number appears in outsite
siteRow <- exBySite[[i]]
if(i %in% outSite){
##extracts the data from the site-pair table by site
siteRow <- siteRow[grep("s2\\.", colnames(siteRow))]
colnames(siteRow) <- sapply(strsplit(colnames(siteRow), "s2\\."), "[[", 2)
}else{
##extracts the data from the site-pair table by site
siteRow <- siteRow[grep("s1\\.", colnames(siteRow))]
colnames(siteRow) <- sapply(strsplit(colnames(siteRow), "s1\\."), "[[", 2)
}
exBySite[[i]] <- siteRow
}
##transforms data from list to table
siteData <- do.call("rbind", exBySite)
##sets up objects to be returned by the function
modelTestValues <- matrix(NA,4,nVars*5,dimnames = list(c("Model deviance",
"Percent deviance explained",
"Model p-value",
"Fitted permutations"),
c("All predictors",
"1-removed",
paste(seq(2,nVars*5-1), "-removed", sep=""))))
##deviance reduction predictor table
varImpTable <- matrix(NA, nVars, nVars*5-1)
rownames(varImpTable) <- varNames
colnames(varImpTable) <- c("All predictors",
"1-removed",
paste(seq(2,nVars*5-2), "-removed", sep=""))
##p value predictor table
pValues <- nModsConverge <- varImpTable
##assigns given site-pair table to new variable, to prevent changing the original input
currSitePair <- spTable
nullGDMFullFit <- 0 ##a variable to track if the fully fitted gdm model returned a NULL object
# create the set up permuted site-pair tables to be used for all
# downstream analyses
message(paste0("Creating ", nPerm, " permuted site-pair tables..."))
# use replicate if number of perms is small or parallel == FALSE
if(parallel==F | nPerm <= 25){
permSpt <- pbreplicate(nPerm, list(permutateSitePair(currSitePair,
siteData,
indexTab,
varNames)))}
# create permuted tables in parallel otherwise
if(parallel==T & nPerm > 25){
# set up parallel processing
cl <- makeCluster(cores)
registerDoParallel(cl)
permSpt <- foreach(k=1:nPerm,
.verbose=F,
.packages=c("gdm"),
.export=c("permutateSitePair")) %dopar%
permutateSitePair(currSitePair, siteData, indexTab, varNames)
stopCluster(cl)}
# create new vector to track varNames
varNames.x <- varNames
message("Starting model assessment...")
for(v in 1:length(varNames)){
# ends the loop if only 1 variable remains
if(length(varNames.x)<2){
stop("Only one predictor remains...variable assessment stopped.")
}
if(v>1){
message(paste0("Removing ", names(elimVar), " and proceeding with the next round of permutations."))
}
if(is.numeric(splines)){
splines <- rep(unique(splines), length(varNames.x))
}
# runs gdm, first time on site-pair table with all variables
# As variables are eliminated the "full" site-pair table will
# contain fewer variables
fullGDM <- gdm(currSitePair, geo=geo, splines=splines, knots=knots)
message(paste0("Percent deviance explained by the full model = ", round(fullGDM$explained,3)))
if(is.null(fullGDM)==TRUE){
warning(paste("The model did not converge when testing variable: ", varNames.x[v],
". Terminating analysis and returning output completed up to this point.", sep=""))
break
}
# fit gdm to each permuted site-pair table
message("Fitting GDMs to the permuted site-pair tables...")
permGDM <- lapply(permSpt, function(x){
gdm(x, geo=geo, splines=splines, knots=knots)
})
##extracts deviance of permuted gdms
permModelDev <- sapply(permGDM, function(mod){mod$gdmdeviance})
##if needed, removes nulls from output permModelDev
modPerms <- length(which(sapply(permModelDev,is.null)==TRUE))
if(modPerms>0){
permModelDev <- unlist(permModelDev[-(which(sapply(permModelDev,is.null)==T))])
}
if(v>1){
colnames(modelTestValues)[v] <- colnames(pValues)[v] <- colnames(varImpTable)[v] <- colnames(nModsConverge)[v] <- paste0(names(elimVar),"-removed")
}
##begins to fill in the output table with data from fully fitted model
modelTestValues[1,v] <- round(fullGDM$gdmdeviance,3)
modelTestValues[2,v] <- round(fullGDM$explained,3)
#p-value
modelTestValues[3,v] <- round(sum(permModelDev<=fullGDM$gdmdeviance)/(nPerm-modPerms),3)
##fitted permutations
modelTestValues[4,v] <- round(nPerm-modPerms,0)
# now permutate each variable individually and fit models
if(parallel==TRUE){
if(length(varNames.x)<cores){
cores <- length(varNames.x)
}
# set up parallel processing
cl <- makeCluster(cores)
registerDoParallel(cl)
# foreach function to create site-pair tables with each variable permuted,
# fit gdms and extract deviance.
permVarDev <- foreach(k=1:length(varNames.x),
.verbose=F,
.packages=c("gdm"),
.export = c("currSitePair")) %dopar%{
if(varNames.x[k]!="Geographic"){
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
idx1 <- grep(paste("^s1\\.", varNames.x[k], "$", sep=""), colnames(x))
idx2 <- grep(paste("^s2\\.", varNames.x[k], "$", sep=""), colnames(x))
spt[,c(idx1, idx2)] <- x[,c(idx1, idx2)]
return(spt)
})}
if(varNames.x[k]=="Geographic"){
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
s1 <- sample(1:nrow(spt), nrow(spt))
s2 <- sample(1:nrow(spt), nrow(spt))
s3 <- sample(1:nrow(spt), nrow(spt))
s4 <- sample(1:nrow(spt), nrow(spt))
spt[,3] <- spt[s1,3]
spt[,4] <- spt[s2,4]
spt[,5] <- spt[s3,5]
spt[,6] <- spt[s4,6]
return(spt)
})}
gdmPermVar <- lapply(lll, function(x){
try(gdm(x, geo=geo, splines=splines, knots=knots))
})
##extracts deviance of permuted gdms
permModelDev <- sapply(gdmPermVar, function(mod){mod$gdmdeviance})
return(permModelDev)
}
##closes cores
#close(pb)
stopCluster(cl)
}
if(parallel==FALSE){
# for-loop to create site-pair tables with each variable permuted,
# fit gdms and extract deviance.
permVarDev <- list()
for(k in 1:length(varNames.x)){
if(varNames.x[k]!="Geographic"){
message(paste0("Assessing importance of ", varNames.x[k], "..."))
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
idx1 <- grep(paste("^s1\\.", varNames.x[k], "$", sep=""), colnames(x))
idx2 <- grep(paste("^s2\\.", varNames.x[k], "$", sep=""), colnames(x))
spt[,c(idx1, idx2)] <- x[,c(idx1, idx2)]
return(spt)
})}
if(varNames.x[k]=="Geographic"){
message("Assessing importance of geographic distance...")
# permute a single variable
lll <- lapply(permSpt, function(x, spt=currSitePair){
s1 <- sample(1:nrow(spt), nrow(spt))
s2 <- sample(1:nrow(spt), nrow(spt))
s3 <- sample(1:nrow(spt), nrow(spt))
s4 <- sample(1:nrow(spt), nrow(spt))
spt[,3] <- spt[s1,3]
spt[,4] <- spt[s2,4]
spt[,5] <- spt[s3,5]
spt[,6] <- spt[s4,6]
return(spt)
})}
gdmPermVar <- lapply(lll, function(x){
try(gdm(x, geo=geo, splines=splines, knots=knots))
})
##extracts deviance of permuted gdms
permVarDev[[k]] <- unlist(sapply(gdmPermVar, function(mod){mod$gdmdeviance}))
}
}
names(permVarDev) <- varNames.x
nullDev <- fullGDM$nulldeviance
for(var in varNames.x){
grepper <- grep(paste0("^",var,"$"), names(permVarDev))
varDevTab <- permVarDev[[grepper]]
# number of perms for which GDM converged
nConv <- length(varDevTab)
nModsConverge[which(rownames(varImpTable) == var),v] <- nConv
# remove NULLs (GDMs that did not converge)
#if(nConv>0){
#varDevTab <- unlist(varDevTab[-(which(sapply(is.null(varDevTab))))])
#}
# calculate variable importance
varDevExplained <- 100*(nullDev-varDevTab)/nullDev
varImpTable[which(rownames(varImpTable) == var),v] <- median(100 * abs((varDevExplained - fullGDM$explained)/fullGDM$explained))
if(var!="Geographic"){
##select variable columns to be removed from original site-pair table
testVarCols1 <- grep(paste("^s1\\.", var, "$", sep=""), colnames(currSitePair))
testVarCols2 <- grep(paste("^s2\\.", var, "$", sep=""), colnames(currSitePair))
testSitePair <- currSitePair[,-c(testVarCols1, testVarCols2)]
##run gdm for the missing variable
noVarGDM <- gdm(testSitePair, geo=geo, splines=splines[-1], knots=knots)
} else {
noVarGDM <- gdm(currSitePair, geo=F, splines=splines[-1], knots=knots)
}
# calculate p-Value
permDevReduct <- noVarGDM$gdmdeviance - varDevTab
pValues[which(rownames(pValues) == var),v] <- sum(permDevReduct>=(varDevTab - fullGDM$gdmdeviance))/(nConv)
}
if(max(na.omit(pValues[,v]))<pValue){
message("All remaining predictors are significant, ceasing assessment.")
message(paste0("Percent deviance explained by final model = ", round(fullGDM$explained,3)))
message("Final set of predictors returned: ")
for(vvv in 1:length(fullGDM$predictors)){
message(fullGDM$predictors[vvv])
}
break
}
if(predSelect==T){
# eliminate least important variable
#elimVar <- which.min(varImpTable[,v])
# eliminate variable with highest p-value
elimVar <- which.max(pValues[,v])
if(names(elimVar)!="Geographic"){
##select variable columns to be removed from original site-pair table
remVarCols1 <- grep(paste("^s1\\.", names(elimVar), "$", sep=""), colnames(currSitePair))
remVarCols2 <- grep(paste("^s2\\.", names(elimVar), "$", sep=""), colnames(currSitePair))
# remove columns from site-pair table
currSitePair <- currSitePair[,-c(remVarCols1, remVarCols2)]
# remove columns from permuted site-pair tables
permSpt <- lapply(permSpt, function(x){
x[,-c(remVarCols1, remVarCols2)]
})
} else {
geo <- F
}
varNames.x <- varNames.x[-which(varNames.x==names(elimVar))]
}
if(v==1 & predSelect==F){
message("Backwards elimination not selected by user (predSelect=F). Ceasing assessment.")
message(paste0("Percent deviance explained by final model = ", round(fullGDM$explained,3)))
message("Final set of predictors returned: ")
for(vvv in 1:length(fullGDM$predictors)){
message(fullGDM$predictors[vvv])
}
break
}
}
if(v==1 & predSelect==F){
# Model assessment
modelTestVals <- data.frame(matrix(round(modelTestValues[,1], 3), ncol=1))
rownames(modelTestVals) <- rownames(modelTestValues)
colnames(modelTestVals) <- "All predictors"
#Variable importance
varImpTab <- data.frame(matrix(round(varImpTable[,1], 3), ncol=1))
rownames(varImpTab) <- rownames(varImpTable)
colnames(varImpTab) <- "All predictors"
#Variable selection p-values
pVals <- varImpTab
pVals[,1] <- round(pValues[,1], 3)
#Variable selection model convergence
nModsConv <- varImpTab
nModsConv[,1] <- round(nModsConverge[,1], 3)
outObject <- list(modelTestVals, varImpTab, pVals, nModsConv)
names(outObject) <- c("Model assessment", "Predictor Importance", "Predictor p-values", "Model Convergence")
} else{
outObject <- list(round(modelTestValues[,1:v], 3), round(varImpTable[,1:v],3), round(pValues[,1:v],3), nModsConverge[,1:v])
names(outObject) <- c("Model assessment", "Predictor Importance", "Predictor p-values", "Model Convergence")
}
##if given, writes out files to space on disk
if(is.null(outFile)==FALSE){
save(outObject, file=outFile)
}
return(outObject)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.