R/glmm_wrap.R

In drmarcogir/Rmacdiv:

glmm_wrap<-function(fixed,random,indf,scale,group,method,dist,sitecol){
    options(na.action = "na.fail")
    # relative importance: Akaike weights
    importance<-NULL
    coefficients<-NULL
    residuals<-NULL
    # residuals with model ranking
    scale.l<-unique(indf[,scale])
    # loop through scale
    for (i in 1:length(scale.l)){
        tmp<-subset(indf,scale==scale.l[i])
        sp.group.l<-as.character(unique(tmp[,group]))
        # loop through groups
        for (j in 1:length(sp.group.l)){
            tmp1<-subset(tmp,sp.group==sp.group.l[j])
            # fit full model
            inform<-as.formula(paste(fixed,"+(1|",random,")",sep=""))
            if(dist!="gaussian"){
                mod<-glmer(inform,data=tmp1,family=dist)
            } else {
                mod<-lmer(inform,data=tmp1,REML=FALSE)
            }
            # model selection
            modt<-dredge(mod)
            # model averaging
            modavgest <- model.avg(modt, cumsum(weight) <= .95)
            # get models
            confset <- get.models(modt, subset = cumsum(weight) <= .95)
            # store parameter estimates
            pardf<-data.frame(variable=names(modavgest$coefficients[1,][-1]),parest=modavgest$coefficients[1,][-1],
                              SE=summary(modavgest)$coefmat.full[-1,3],P=summary(modavgest)$coefmat.full[-1,5])
            pardf$scale<-scale.l[i]
            pardf$sp.group<-sp.group.l[j]
            coefficients<-rbind(pardf,coefficients)
            # residuals for best performing model
            resdf<-data.frame(residuals=residuals(confset[row.names(modt)[1]][[1]]),
                              tmp1[,c("Island",sitecol,"scale","sp.group")])
            residuals<-rbind(resdf,residuals)
            if(method=="mumin"){
                rel.imp<-data.frame(variable=names(summary(modavgest)$importance),importance=summary(modavgest)$importance,length=dim(modavgest$msTable)[1])
                rel.imp$scale<-scale.l[i]
                rel.imp$sp.group<-sp.group.l[j]
                importance<-rbind(rel.imp,importance)
            }
            if(method=="cade"){
                # ------- Model averaging: Cade's method (2015)
                # subset only models within confidence set
                modt95<-modt[names(confset),]
                sum.wt<-sum(modt95$weight)
                # get model table with select columns
                as.data.frame(modt95) %>%
                    tibble::rownames_to_column() %>%
                    rename(Int=`(Intercept)`) %>%
                    dplyr::select(-c(Int,df,logLik,AICc,delta)) ->modtb1
                # results dataframe (stores weights and ratios for variables)
                results<-NULL
                # goes through each variable in turn
                for (z in 2:c(dim(modtb1)[2]-1)){
                    modtb1[,c("rowname",names(modtb1[z]),"weight")] %>%
                        rename("var" = !!names(.[2])) %>%
                        dplyr::filter(!is.na(var)) ->varint
                    confset[varint$rowname] ->modlist
                    for(a in 1:length(modlist)){
                        # extract coefficients
                        coefdf<-data.frame(variable=row.names(coef(summary(modlist[[a]]))),coef(summary(modlist[[a]])))
                        row.names(coefdf)<-1:dim(coefdf)[1]
                        names(coefdf)[4]<-c("z.value")
                        # get t-score for variable of interest
                        coefdf %>%
                            dplyr::select(c(variable,z.value)) %>%
                            filter(variable==names(modtb1[z]))  %>%
                            pull(z.value) ->t.varint
                        # absolute value for variable of interest
                        t.varint<-abs(t.varint)
                        # get t-scores for other variables
                        dplyr::select(coefdf,c(variable,z.value)) %>%
                            filter(variable!=names(modtb1[z]) & variable!="(Intercept)") %>%
                            pull(z.value) ->t.other
                        if(length(t.other)==0){
                            finalres<-data.frame(variable=names(modtb1[z]),ratio=t.varint/t.varint,
                                                 weight=filter(modtb1,rowname==names(modlist)[a])$weight)
                        } else {
                            # max absolute value for other variables
                            t.other<-unique(max(abs(t.other)))
                            # final dataframe with model weights
                            finalres<-data.frame(variable=names(modtb1[z]),
                                                 ratio=t.varint/t.other,weight=filter(modtb1,rowname==names(modlist)[a])$weight)
                        }
                        # bind results
                        results<-rbind(finalres,results)
                    }
                }
                # get variable importance
                results %>%
                    mutate(rmw=ratio*as.numeric(weight)) %>%
                    dplyr::select(c(variable,rmw)) %>%
                    group_by(variable) %>%
                    summarise(importance=sum(rmw)) %>%
                    mutate(importance=importance/sum.wt) %>%
                    as.data.frame() ->rel.imp
                # insert information on scale & sp.group
                #rel.imp<-data.frame(variable=names(summary(modavgest)$importance),
                #                   importance=rel.imp)
                rel.imp$scale<-scale.l[i]
                rel.imp$sp.group<-sp.group.l[j]
                importance<-rbind(rel.imp,importance)
            } # end of if statement
        }      # end of j loop
    } # end of i loop
    final.results<-list(coefficients,importance,residuals)
    names(final.results)<-c("coefficients","importance","residuals")
    return(final.results)
}