#R-script for non-linear multi quantile regression and lowest AIC selection require(splines) require(quantreg) #grzq function########################################################################## grzq<-function(dvar,indvar,tau=seq(0.01,0.99,0.005),method="lasso",lambda=1){ # dvar(vector): values of dependent variable # indvar(vector): values of independent variable ###quantreg!!!### #function for quantile regression rqfun<-function(z) rq(log1p(dvar)~bs(indvar,degree=z),tau=tau,method=method,lambda=lambda) ##spline selection #range of spline degrees bs<-c(2:5) #models for different bs degrees rqlist<-lapply(bs,rqfun) #maximal value forcasted for each quantile obsv<-lapply(rqlist,function(z) apply(expm1(fitted(z)),2,max)) #selection of the quantiles foreccasting abundance > than the minimal observed obsv2<-do.call("cbind",lapply(obsv,function(z)(z>=min(dvar[dvar>0])))) #AICc for each quantile AICc<-function(z){ zz<-dim(z[[1]])[1]-1 aicc<-AIC(z)+2*zz*(zz+1)/(length(indvar)-zz-1) return(aicc) } bbs<-sapply(rqlist,AICc) bbs2<-array(c(bbs,obsv2),dim=c(dim(bbs),2)) ##mean AIC for quantiles predicting occurrence bbs3<-apply(bbs2,2,function(z)sum(z[,1]*z[,2])/length(which(z[,2]==T)))) #model with the lowest average AIC rqmod<-rqlist[[which.min(bbs3)]] #bs degree giving the lowest average AIC sdg<-which.min(bbs3)+1 return(list(rqmod,sdg)) \end{document} }