# R/multiSAR.R In mmSAR: mmSAR is an R package for the modelling of the Species-Area Relationship (SAR).

#### Documented in multiSAR

```multiSAR <-
function(modelList,data,nBoot=999,crit="Info",norTest="lillie",alpha=0.05,verb=FALSE) {

######## multiSAR : model selection and averaging function
#
# modelList a vector of model names : ex : c("power","expo","negexpo","monod","logist","ratio","lomolino","weibull")
# data : an mmSAR data object : a list (run data(data.galap)) for an example
# nBoot : the number of bootstrap resamples
# crit : "Info" or "Bayes"
# norTest : "lillie" or "shapiro" or "null"
# alpha : level of tests for normality and homoscedasticity
# verb : FALSE or TRUE
#
##########################################################

#will check if the data is OK
nPoints <- length(data\$data[[1]])
nlig <- length(modelList)

############Criteria must be "Info" for AIC/AICc or "Bayes" for BIC
#choosing an IC criterion (AIC or AICc or BIC)
if(crit == "Info") {
if ( (nPoints / 3) < 40 ) { IC = "AICc" } else { IC = "AIC"}
} else {
if(crit == "Bayes") { IC = "BIC" } else { stop("Criteria must be 'Info' for AIC/AICc or 'Bayes' for BIC")}
}

library(numDeriv)

#GO
if(verb) cat("##### multiSAR #####\n#")
if(verb) cat("# Choosen criterion is ",IC,"\n")

#This should eventually be an option
#Test on data points (if one richness == 0 then the data point is deleted)
#isNull = which(data\$data[[2]]==0)

#if (length(isNull)!=0) {
#	if(verb) cat("Dataset contained ",length(isNull)," zero abundance point(s) that was(were) deleted for analysis\n")
#	data\$data = data\$data[-isNull,]
#}#end of if isNull

#matrix of optimisation results
vars <- c("p1","p2","p3","AICc","D.AICc","AICcW","AIC","D.AIC","AICW","BIC","D.BIC","BICW","RSS","R2","Norm Stat","Norm p.val","Pearson","Pea p.val")
optimResult = matrix(0,nlig,length(vars))
colnames(optimResult) = vars
rownames(optimResult) <- modelList

#List of Jacobian and Hat matrix
matList = list()

#matrix of calculated values and residuals and transformed residuals
pointsNames <- paste("S",c(1:nPoints))
calculated <- residuals <- transResiduals <- matrix(0,nlig,length(pointsNames))
colnames(calculated) <- colnames(residuals) <- colnames(transResiduals) <- pointsNames
rownames(calculated) <- rownames(residuals) <- rownames(transResiduals) <- modelList

#vector of final values (model averaging)
finalVect <- vector("numeric",length(pointsNames))
names(finalVect) <- pointsNames

#getting optimization results for all models
for (i in 1:nlig){

for (j in 1:eval(parse(text=as.character(modelList[i])))\$paramnumber) {optimResult[i,paste("p",j,sep="")] <- optimres\$par[j]}
optimResult[i,"AIC"] <- optimres\$AIC
optimResult[i,"AICc"] <- optimres\$AICc
optimResult[i,"BIC"] <- optimres\$BIC
optimResult[i,"R2"] <- optimres\$R2
optimResult[i,"Norm Stat"] <- optimres\$normaStat
optimResult[i,"Norm p.val"] <- optimres\$normaPval
optimResult[i,"Pearson"] <- optimres\$pearson
optimResult[i,"Pea p.val"] <- optimres\$pearpval

calculated[i,] <- optimres\$calculated
residuals[i,] <- optimres\$calculated - data\$data[,2]

#jacobian and Hat Matrix
#first data Point

for (k in 2:nPoints) {
}

jacobbis <- t(jacob)%*%jacob
s <- svd(jacobbis)
jacobbismun = s\$v%*%(diag(1/s\$d))%*%(t(s\$u))
hatMat = jacob%*%jacobbismun%*%t(jacob)
matList[[i]] <- list(jacob=jacob,hatMat=hatMat)

#Residuals transformation from Davidson and Hinkley, 1997 "Bootstrap methods and their applications" p 259 eq (6.9)
diagHatMat = diag(hatMat)
transResiduals[i,] <- residuals[i,] - mean(residuals[i,])
transResiduals[i,] <- transResiduals[i,] / sqrt( 1 - diagHatMat )

}#end of for

names(matList) = modelList

#Fit validation
flags <- vector("numeric",nlig)

if(norTest!="null"){
for (i in 1:nlig) { if (optimResult[i,"Norm p.val"]<alpha || optimResult[i,"Pea p.val"]<alpha) {flags[i]<-"KO"} else {flags[i]<-"OK"}  }
}else{
flags <- rep("OK",nlig)
}#eo ifelse

filtOptimResult <- subset(optimResult,flags=="OK")
filtCalculated <- subset(calculated,flags=="OK")
filtMatList <- matList[flags=="OK"]
filtModelList <- modelList[flags=="OK"]

if(verb) cat("# Valid models : ", paste(rownames(filtOptimResult),collapse=", "),"\n#\n")

#Models comparison

DeltaICvect <- vector()
akaikeweightvect <- vector()

filtNlig <- dim(filtOptimResult)[1]

for (i in 1:filtNlig){

#Delta IC = ICi - ICmin
DeltaIC <- filtOptimResult[i,IC] - min(filtOptimResult[,IC])
DeltaICvect <- c(DeltaICvect,DeltaIC)
}

for (i in 1:filtNlig){
#Akaike Weigths
akaikesum <- sum(exp( -0.5*(DeltaICvect)))
akaikeweight <- exp(-0.5*DeltaICvect[i]) / akaikesum
akaikeweightvect <- c(akaikeweightvect,akaikeweight)
}

columnDelta = paste("D.",IC,sep="")
filtOptimResult[,columnDelta] <- DeltaICvect
columnW = paste(IC,"W",sep="")
filtOptimResult[,columnW] <- akaikeweightvect

if(verb) cat("# Akaike weigths : ", paste(round(filtOptimResult[,columnW],4),collapse=", "),"\n#\n")

#Averaging
for (i in 1:nPoints) {
finalVect[i] <- sum(akaikeweightvect*filtCalculated[,i])
}

#Averaging validation
avResiduals = data\$data[[2]] - finalVect
shapRes= shapiro.test(avResiduals)
if(verb) cat("# Averaging residuals normality (p.value) : ",shapRes\$p.value,"\n#\n")

cor <- cor.test(avResiduals,data\$data[[1]])
if(verb) cat("# Averaging residuals/X values correlation (method: ",cor\$method,") (Value,p.value) : ",cor\$estimate,",",cor\$p.value,"\n#\n")

################################################################################
#Bootstrapping residuals and model averaging
################################################################################

#Matrix of boot Samples
bootMatrix=matrix(0,nBoot,nPoints)

#array of optimisation results
optimBootResult = array(0,c(nlig,length(vars),nBoot),dimnames=list(modelList,vars,seq(1,nBoot)))

#array of calculated values
bootCalculated <- array(0,c(nlig,length(pointsNames),nBoot),dimnames=list(modelList,pointsNames,seq(1,nBoot)))

#flags for fitting validation
flags <- matrix(0,nlig,nBoot)

if(verb) cat("# Bootstrap Samples creation\n#\n")

#vector of choosen models
choosenModels = vector()

#test variable
nGoodBoot = 1

while (nGoodBoot < nBoot+1) {

test <- 1

chousModel = filtModelList[rmultinom(1, 1, akaikeweightvect)==1]
if(verb) cat("Boot Sample : ",nGoodBoot," Choose model : ",chousModel,"\n")

choosenModels[nGoodBoot] = chousModel

while (test != 0 ) {

for (l in 1:nPoints) {
positives = transResiduals[chousModel,][transResiduals[chousModel,] > 0]
negatives = transResiduals[chousModel,][transResiduals[chousModel,] < 0]

if (calculated[chousModel,l] > 0 ) {
vtci = negatives[abs(negatives) <= calculated[chousModel,l] ]
vtci = c(vtci,positives)
value = sample(vtci,1)
} else {

vtci = positives[positives >= abs(calculated[chousModel,l]) ]
value = sample(vtci,1)
}

bootMatrix[nGoodBoot,l] <- calculated[chousModel,l] + value
}#end of for

#test if one species richness is negative
test=length( which(bootMatrix[nGoodBoot,]<0) )

if(verb) cat("BootSample : ",bootMatrix[nGoodBoot,],"\n")

}#end of while

###########################################################
#Do the model averaging for each bootstrap sample
###########################################################

if(verb) cat("# model optimization on bootstrap resamples\n#\n")

for (k in 1:nlig){

optimres = tryCatch(rssoptim(eval(parse(text=as.character(modelList[k]))),data=list(name="bootSample",data=data.frame(a=data\$data[[1]],s=bootMatrix[nGoodBoot,])),norTest,verb),error = function(e) { list(convergence=999) } )

if (optimres\$convergence != 0) {
} else {
for (j in 1:eval(parse(text=as.character(modelList[k])))\$paramnumber){optimBootResult[k,paste("p",j,sep=""),nGoodBoot] <- optimres\$par[j]}
optimBootResult[k,"AIC",nGoodBoot] <- optimres\$AIC
optimBootResult[k,"AICc",nGoodBoot] <- optimres\$AICc
optimBootResult[k,"BIC",nGoodBoot] <- optimres\$BIC
optimBootResult[k,"R2",nGoodBoot] <- optimres\$R2
optimBootResult[k,"Norm Stat",nGoodBoot] <- optimres\$normaStat
optimBootResult[k,"Norm p.val",nGoodBoot] <- optimres\$normaPval
optimBootResult[k,"Pearson",nGoodBoot] <- optimres\$pearson
optimBootResult[k,"Pea p.val",nGoodBoot] <- optimres\$pearpval
bootCalculated[k,,nGoodBoot] <- optimres\$calculated

#Fitting validation
if(norTest!="null"){
if (optimBootResult[k,"Norm p.val",nGoodBoot]<alpha || optimBootResult[k,"Pea p.val",nGoodBoot]<alpha || length(which(bootCalculated[k,,nGoodBoot]<0)) !=0 ) { flags[k,nGoodBoot]<-"KO"
} else {
flags[k,nGoodBoot]<-"OK"
}#end of if/else on Shap and Corr
}else{
flags[k,nGoodBoot]<-"OK"
}#eo if/else norTest!=NULL
}#end of if/else on convergence 2
}#end of if/else on convergence 1

}#end of for k

if ( length(which(flags[,nGoodBoot]!="KO")) == 0 ) badBoot=TRUE
if ( length(which(flags[,nGoodBoot]==0)) != 0 ) badBoot=TRUE

#write the bootSample to a file
#bootFileName = paste("bootSamp_",data\$name,".txt",sep="")
#bootText = paste("BootSample",nGoodBoot,"\n",sep="")
#cat(bootText,file = bootFileName,append=TRUE)
#write(bootMatrix[nGoodBoot,], file = bootFileName,ncolumns= nPoints, append = TRUE, sep = " ")
nGoodBoot <- nGoodBoot + 1
}#end of while

#Applying the filter (flags)
#transform 3D table to list
filtOptimBootResult=vector("list", nBoot)
for (i in 1:nBoot) filtOptimBootResult[[i]] <- optimBootResult[,,i]
for (i in 1:nBoot) filtOptimBootResult[[i]] <- subset(filtOptimBootResult[[i]],flags[,i]=="OK")
filtBootCalculated = vector("list", nBoot)
for (i in 1:nBoot) filtBootCalculated[[i]] <- bootCalculated[,,i]
for (i in 1:nBoot) filtBootCalculated[[i]] <- subset(filtBootCalculated[[i]],flags[,i]=="OK")

bootHat = matrix(0,nBoot,nPoints)
f=0

if(verb) cat("# model averaging on bootstrap resamples\n#\n")

for (k in 1:nBoot) {

DeltaICvect <- vector()
akaikeweightvect <- vector()
filtNlig <- dim(filtOptimBootResult[[k]])[1]

if (filtNlig != 0) {

for (i in 1:filtNlig){

#Delta IC = ICi - ICmin
DeltaIC <- filtOptimBootResult[[k]][i,IC] - min(filtOptimBootResult[[k]][,IC])
DeltaICvect <- c(DeltaICvect,DeltaIC)
}#end of for i

for (i in 1:filtNlig){
#Akaike Weigths
akaikesum <- sum(exp( -0.5*(DeltaICvect)))
akaikeweight <- exp(-0.5*DeltaICvect[i]) / akaikesum
akaikeweightvect <- c(akaikeweightvect,akaikeweight)
}#end of for i

columnDelta = paste("D.",IC,sep="")
columnW = paste(IC,"W",sep="")

filtOptimBootResult[[k]][,columnDelta] <- DeltaICvect
filtOptimBootResult[[k]][,columnW] <- akaikeweightvect

#Averaging
for (i in 1:nPoints) {
bootHat[k,i] <- sum(akaikeweightvect*filtBootCalculated[[k]][,i])
}#end of for i

} else { bootHat[k,] <- rep(0,nPoints) }

} #end of for k

bootSort=apply(bootHat,2,sort)

res=list(data=data,models=modelList,optimRes=optimResult,filtOptimRes=filtOptimResult,calculated=calculated,filtCalculated=filtCalculated,averaged=finalVect,DeltaIC=DeltaICvect,akaikeweight=akaikeweightvect,avResiduals=avResiduals,shapAvRes=shapRes,corAvRes=cor,bootMatrix=bootMatrix,optimBootResult=optimBootResult,bootCalculated=bootCalculated,flags=flags,filtOptimBootResult=filtOptimBootResult,filtBootCalculated=filtBootCalculated,bootSort=bootSort,bootHat=bootHat,bootMatrix=bootMatrix,IC=IC)

if(verb) cat("###################\n")

invisible(res)

} #end of multiSAR
```

## Try the mmSAR package in your browser

Any scripts or data that you put into this service are public.

mmSAR documentation built on May 31, 2017, 3:27 a.m.