Rmacdiv:

#' Function for using a series of mantel tests using multivariate matrices
#' methods of Carvalho et al. (2012)
#'
#' @ insp = matrix containing species community data
#' @ inenv = list containing environmental data
#' @ group =character vector indicating column name for site IDs (if there is more than one site )
#'
#'


mrmfit<-function(spdat,topodat=NULL,lusedat=NULL){
    # ---first step: spider dissimilarity matrices
    spdist<-vector("list")
    # list to loop through
    comp.l<-unique(spdat$comp)
    isl.l<-unique(spdat$site)
    hab.l<-unique(spdat$habitat)
    for (i in 1:3){
        # subset by biodiversity component
        t1<-subset(decayres,comp==comp.l[i])
        for (j in 1:length(isl.l)){
            # subset by island
            t2<-subset(t1,site==isl.l[j])
            for (z in 1:length(hab.l)){
                # subset by habitat
                t3<-subset(t2,habitat==hab.l[z])
                # new plot names
                t3$p1<-stri_split_fixed(t3$comparison,"-",simplify=T)[,1]
                t3$p2<-stri_split_fixed(t3$comparison,"-",simplify=T)[,2]
                # dissimilarity matrix
                betavals <-t3[,c("beta")]
                nams <- with(t3, unique(c(as.character(p1), as.character(p2))))
                attributes(betavals) <- with(t3, list(Size = length(nams),
                                                      Labels = nams,
                                                      Diag = FALSE,
                                                      Upper = FALSE,
                                                      method = "user"))
                class(betavals) <- "dist"
                betavals1<-list(betavals)
                names(betavals1)<-paste(comp.l[i],"_",isl.l[j],"_",gsub(" ",".",hab.l[z]),sep="")
                spdist<-append(spdist,betavals1)
            }
        }
    }

    # --- second step: topographic dissimilarity between sites

    topodist<-vector("list")

    for (i in 1:length(topodat)){
        # subset scale
        topotmp<-topodat[[i]]
        # subset island
        isl.l<-as.character(unique(topotmp$Island))
        for (j in 1:length(isl.l)){
            topotmp1<-subset(topotmp,Island==isl.l[j])
            # assign new row number
            row.names(topotmp1)<-(topotmp1$Database.code)
            # create distance matrix
            topotmp2<-vegdist(topotmp1[,c("alt_range","alt_sd","aspect_sd",
                                          "slope_range","slope_sd")],method="jaccard")
            topotmp3<-list(topotmp2)
            names(topotmp3)<-paste(isl.l[j],"_",names(topodat)[i],sep="")
            topodist<-append(topodist,topotmp3)
        }

    }

    # --- third step: land use dissimilarity between sites
    if(is.null(lusedat)){
    } else {
      lusedist<-vector("list")
      for (i in 1:length(lusedat)){
        # subset scale
        lusetmp<-lusedat[[i]]
        # loop through islands
        for (j in 1:length(lusetmp)){
          lusetmp1<-select(lusetmp[[j]],-c(Island,Database.code))
          lusetmp2<-vegdist(lusetmp1,method="jaccard")
          lusetmp3<-list(lusetmp2)
          names(lusetmp3)<-paste(names(lusetmp)[j],"_",names(lusedat)[i],sep="")
          lusedist<-append(lusedist,lusetmp3)
        }
      }
    }

    # --- fourth step: MRM models
    # topo + luse models

    if(!is.null(lusedat)){

        modelres<-NULL

        for (i in 1:length(spdist)){
            # go through each island
            # split names
            mydf<-stri_split_fixed(names(spdist)[i],"_",simplify=TRUE)
            # match island to topography and land use data
            toponames<-names(topodist)[str_detect(names(topodist),mydf[,2])]
            lusenames<-names(lusedist)[str_detect(names(lusedist),mydf[,2])]
            # extract data
            # topography info
            topodist.isl<-topodist[toponames]
            if(length(lusenames)==0){
            } else {
                # land use info
                lusedist.isl<-lusedist[lusenames]
            }
            # fit models go through each scale
            for (j in 1:6){
                # uses topography as an index (for scale!)
                myind<-stri_split_fixed(names(topodist.isl)[j],"_",simplify=TRUE)[,2]
                # predictors
                scaleinfo<-as.data.frame(stri_split_fixed(names(topodist.isl),"_",simplify=TRUE))
                scaleinfo1<-scaleinfo[scaleinfo$V2 %in% myind,]
                predname<-paste(scaleinfo1$V1,"_",scaleinfo1$V2,sep="")
                topo<-topodist.isl[[predname]]
                luse<-lusedist.isl[[predname]]
                # fit model
                if(is.null(luse)){
                    mod<-try(MRM(spdist[[i]]~topo, nperm=30000),silent=TRUE)
                } else {
                    mod<-try(MRM(spdist[[i]]~topo+luse, nperm=30000),silent=TRUE)
                }

                if(class(mod)=="try-error"){
                    dfres<-data.frame(variable=c("topo","luse"),coef=NA,p.coef=NA,R2=NA,R2.pvalue=NA,
                                      F.value=NA,F.pval=NA,scale=scaleinfo1$V2,Island=scaleinfo1$V1,habitat=mydf[,3],
                                      comp=mydf[,1])
                    modelres<-rbind(dfres,modelres)

                } else {
                    dfres<-data.frame(variable=names(mod$coef[,1][-1]),coef=mod$coef[,1][-1],p.coef=mod$coef[,2][-1],
                                      R2=mod$r.squared[1],R2.pvalue=mod$r.squared[2],F.value=mod$F.test[1],
                                      F.pval=mod$F.test[2],scale=scaleinfo1$V2,Island=scaleinfo1$V1,habitat=mydf[,3],
                                      comp=mydf[,1])
                    modelres<-rbind(dfres,modelres)
                }
            }

        }
    }

    # topo models only
    if(is.null(lusedat)){

        modelres<-NULL

        for (i in 1:length(spdist)){
            # go through each island
            # split names
            mydf<-stri_split_fixed(names(spdist)[i],"_",simplify=TRUE)
            # match island to topography and land use data
            toponames<-names(topodist)[str_detect(names(topodist),mydf[,2])]
            # extract data
            # topography info
            topodist.isl<-topodist[toponames]
            # fit models go through each scale
            for (j in 1:6){
                # uses topography as an index (for scale!)
                myind<-stri_split_fixed(names(topodist.isl)[j],"_",simplify=TRUE)[,2]
                # predictors
                scaleinfo<-as.data.frame(stri_split_fixed(names(topodist.isl),"_",simplify=TRUE))
                scaleinfo1<-scaleinfo[scaleinfo$V2 %in% myind,]
                predname<-paste(scaleinfo1$V1,"_",scaleinfo1$V2,sep="")
                topo<-topodist.isl[[predname]]
                # fit model
                mod<-try(MRM(spdist[[i]]~topo, nperm=30000),silent=TRUE)
                if(class(mod)=="try-error"){
                  dfres<-data.frame(variable=c("topo"),coef=NA,p.coef=NA,R2=NA,R2.pvalue=NA,
                                    F.value=NA,F.pval=NA,scale=scaleinfo1$V2,Island=scaleinfo1$V1,habitat=mydf[,3],
                                    comp=mydf[,1])
                    modelres<-rbind(dfres,modelres)

                } else {
                  dfres<-data.frame(variable=names(mod$coef[,1][-1]),coef=mod$coef[,1][-1],p.coef=mod$coef[,2][-1],
                                    R2=mod$r.squared[1],R2.pvalue=mod$r.squared[2],F.value=mod$F.test[1],
                                    F.pval=mod$F.test[2],scale=scaleinfo1$V2,Island=scaleinfo1$V1,habitat=mydf[,3],
                                    comp=mydf[,1])
                    modelres<-rbind(dfres,modelres)
                }
            }

        }
    }

    return(modelres)
}