new_multifunc.R

library(multifunc)
data(all_biodepth)
allVars<-qw(biomassY3, root3, N.g.m2,  light3, N.Soil, wood3, cotton3)

germany<-subset(all_biodepth, all_biodepth$location=="Germany")

vars<-whichVars(germany, allVars)

#re-normalize N.Soil so that everything is on the same sign-scale (e.g. the maximum level of a function is the "best" function)
germany$N.Soil<- -1*germany$N.Soil +max(germany$N.Soil, na.rm=TRUE)

germanyThresh<-getFuncsMaxed(germany, vars, threshmin=0.05, threshmax=0.99, prepend=c("plot","Diversity"), maxN=7)
germanyLinearSlopes<-getCoefTab(funcMaxed ~ Diversity, data=germanyThresh, coefVar="Diversity", family=quasipoisson(link="identity"))
germanyLinearIntercepts<-getCoefTab(funcMaxed ~ Diversity, data=germanyThresh, coefVar="(Intercept)", family=quasipoisson(link="identity"))

###### Illustrate Concept
germanyThresh$percent <- 100*germanyThresh$thresholds
myplot <- ggplot(data=germanyThresh, aes(x=Diversity, y=funcMaxed, group=percent)) +
  ylab(expression("Number of Functions" >= Threshold)) +
  xlab("Species Richness") +
  stat_smooth(method="glm", family=quasipoisson(link="identity"), lwd=0.8, fill=NA, aes(color=percent), alpha=0.3) +
  theme_bw(base_size=18) +
  scale_color_gradientn(name="Percent of \nMaximum", colours=rainbow(length(germanyThresh$percent))) 

myplot


myplot+
  geom_hline(yintercept=5,  color="black", lwd=3, lty=2)

myplot+
  geom_hline(yintercept=5,  color="black", lwd=3, lty=2) +
  geom_hline(yintercept=2, color="black", lwd=3, lty=2) 
  

getGerman <- function(div) germanyLinearSlopes$Estimate*div+germanyLinearIntercepts$Estimate

threshAt <- function(nfunc=5, funcs, threshes){
  idx <- max(which(funcs>=nfunc))
  if(is.infinite(idx)) return(0)
  threshes[idx]
}

maxThreshData <- ddply(data.frame(Diversity=1:max(germany$Diversity)), .(Diversity), function(adf){
  funcs <- getGerman(adf$Diversity)
  threshes <- germanyLinearSlopes$thresholds
  data.frame(nfunc=paste(5:3, "Functions", sep=" "),
             maxThresh=c(threshAt(5, funcs, threshes), threshAt(4, funcs, threshes), threshAt(3, funcs, threshes))
  )
})



library(ggplot2)

qplot(Diversity, maxThresh, data=maxThreshData, facets=~nfunc, geom=c("line", "point")) +
  theme_bw(base_size=17) +
  ylab("Maximum Threshold Achieved by n number of functions\n") +
  xlab("\nSpecies Richness")


maxThreshData[which(maxThreshData$nfunc=="5 Functions"),]