[Mmsar-commits] r9 - / pkg pkg/R pkg/data pkg/man www

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Sep 22 14:32:47 CEST 2009


Author: fguilhaumon
Date: 2009-09-22 14:32:47 +0200 (Tue, 22 Sep 2009)
New Revision: 9

Added:
   pkg/
   pkg/DESCRIPTION
   pkg/R/
   pkg/R/backLink.R
   pkg/R/multiSAR.R
   pkg/R/rssoptim.R
   pkg/R/transLink.R
   pkg/data/
   pkg/data/asymp.rda
   pkg/data/data.arr.rda
   pkg/data/data.atl.rda
   pkg/data/data.galap.rda
   pkg/data/data.glea.rda
   pkg/data/data.gleas.rda
   pkg/data/expo.rda
   pkg/data/logist.rda
   pkg/data/lomolino.rda
   pkg/data/monod.rda
   pkg/data/negexpo.rda
   pkg/data/power.rda
   pkg/data/ratio.rda
   pkg/data/weibull.rda
   pkg/man/
   pkg/man/asymp.Rd
   pkg/man/data.arr.Rd
   pkg/man/data.atl.Rd
   pkg/man/data.galap.Rd
   pkg/man/data.glea.Rd
   pkg/man/data.gleas.Rd
   pkg/man/expo.Rd
   pkg/man/logist.Rd
   pkg/man/lomolino.Rd
   pkg/man/mmSAR.Rd
   pkg/man/monod.Rd
   pkg/man/multiSAR.Rd
   pkg/man/negexpo.Rd
   pkg/man/power.Rd
   pkg/man/ratio.Rd
   pkg/man/rssoptim.Rd
   pkg/man/weibull.Rd
   www/
   www/index.php
Log:
Complete package update due to svn missuse.


Added: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	                        (rev 0)
+++ pkg/DESCRIPTION	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,10 @@
+Package: mmSAR
+Type: Package
+Title: mmSAR is an R package for the modelling of the Species-Area Relationship (SAR).
+Version: 0.1
+Date: 2009-09-01
+Author: François Guilhaumon
+Maintainer: François Guilhaumon <francois.guilhaumon at univ-montp2.fr>
+Description: The mmSAR package implement the multimodel SAR paradigm and provides ecologists with user friendly functions fit SAR models and calculate multimodel SAR inference.
+License: GPL
+LazyLoad: yes

Added: pkg/R/backLink.R
===================================================================
--- pkg/R/backLink.R	                        (rev 0)
+++ pkg/R/backLink.R	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,30 @@
+backLink <-
+function(x,boundType) {
+
+	# fonction reciproque de la fonction logit
+	invlogit=function(x) {
+		1/(1+exp(-x))
+	}#end of invlogit
+
+if (length(x) ==1) {
+	res = switch(boundType,
+		R = x,
+		Rplus = exp(x),
+		unif = invlogit(x)
+	)#end of switch
+} else {
+
+res=vector()
+
+	for (i in 1:length(x)) {
+
+		res[i] = switch(boundType[i],
+				R = x[i],
+				Rplus = exp(x[i]),
+				unif = invlogit(x[i])
+			 )#end of switch
+	}#end of for
+res
+}#end of if/else
+}#end of backLink
+

Added: pkg/R/multiSAR.R
===================================================================
--- pkg/R/multiSAR.R	                        (rev 0)
+++ pkg/R/multiSAR.R	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,332 @@
+multiSAR <-
+function(modelList=c("power","expo","negexpo","monod","logist","ratio","lomolino","weibull"),data=gleason,nBoot=999,verb=T,crit="Bayes") {
+
+############Criteria must be "Info" for AIC/AICc or "Bayes" for BIC
+nlig <- length(modelList)
+
+#Test on data points (if one richness == 0 then the data point is deleted)
+isNull = which(data$data[[2]]==0)
+
+if (length(isNull)!=0) {
+	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","R2a","Norm Stat","Norm p.val","Pearson","Pea p.val") #,"t Test","t p.val"
+dig=6 #digits
+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
+nPoints <- length(data$data[[1]])
+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
+
+
+for (i in 1:nlig){
+
+	optimres = rssoptim(eval(parse(text=as.character(modelList[i]))),data,"lillie",graph,verb)
+
+	for (j in 1:eval(parse(text=as.character(modelList[i])))$paramnumber) {optimResult[i,paste("p",j,sep="")] <- round(optimres$par[j],digits=dig)}
+	optimResult[i,"AIC"] <- round(optimres$AIC,digits=dig)
+	optimResult[i,"AICc"] <- round(optimres$AICc,digits=dig)
+	optimResult[i,"BIC"] <- round(optimres$BIC,digits=dig)
+	optimResult[i,"RSS"] <- round(optimres$value,digits=dig)
+	optimResult[i,"R2a"] <- round(optimres$R2a,digits=dig)
+	optimResult[i,"Norm Stat"] <- round(optimres$normaStat,digits=dig)
+	optimResult[i,"Norm p.val"] <- round(optimres$normaPval,digits=dig)
+	optimResult[i,"Pearson"] <- round(optimres$pearson,digits=dig)
+	optimResult[i,"Pea p.val"] <- round(optimres$pearpval,digits=dig)
+
+	calculated[i,] <- optimres$calculated
+	residuals[i,] <- optimres$calculated - data$data[,2]
+	
+	#jacobian and Hat Matrix
+	
+	#first data Point
+	jacob = jacobian( eval(parse(text=as.character(modelList[i])))$rssfun,optimres$par,data=data$data[1,],opt=F)
+
+	for (k in 2:nPoints) {
+	jacob = rbind(jacob,jacobian( eval(parse(text=as.character(modelList[i])))$rssfun,optimres$par,data=data$data[k,],opt=F))
+	}
+	
+	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
+
+
+#Fitting validation
+flags <- vector("numeric",nlig)
+
+for (i in 1:nlig) { if (optimResult[i,"Norm p.val"]<0.05 || optimResult[i,"Pea p.val"]<0.05) {flags[i]<-"KO"} else {flags[i]<-"OK"}  } #|| optimResult[i,"t p.val"]<0.05
+
+filtOptimResult <- subset(optimResult,flags=="OK")
+filtCalculated <- subset(calculated,flags=="OK")
+filtMatList <- matList[flags=="OK"]
+filtModelList <- modelList[flags=="OK"]
+
+#Models comparaison
+
+#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")}
+	}
+
+if(verb) cat("Choosen criterion is ",IC,"\n")
+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] <- round(DeltaICvect,digits=dig)
+#filtOptimResult[,"AICcW"] <- round(akaikeweightvect,digits=5)
+columnW = paste(IC,"W",sep="")
+filtOptimResult[,columnW] <- akaikeweightvect
+
+#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")
+
+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")
+
+################################################################################
+#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("********************************************\n")
+if(verb) cat(" Bootstrap Samples creation and\n")
+if(verb) cat(" Model averaging on boot samples\n")
+if(verb) cat("********************************************\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
+	###########################################################
+
+	for (k in 1:nlig){
+	###########################################################
+	#ICI le tryCatch
+	###########################################################
+
+	badBoot = F
+
+	optimres = tryCatch(rssoptim(eval(parse(text=as.character(modelList[k]))),data=list(name="bootSample",data=data.frame(a=data$data[[1]],s=bootMatrix[nGoodBoot,])),"lillie",graph,verb),error = function(e) {cat("Error from optim function, Swap the bootSample\n") ; list(convergence=999) } )
+
+	if (optimres$convergence != 0) {
+		badBoot=T
+	} else { 
+
+		if (sum(optimres$calculated)==0) { badBoot=T } else { 
+			for (j in 1:eval(parse(text=as.character(modelList[k])))$paramnumber){optimBootResult[k,paste("p",j,sep=""),nGoodBoot] <- round(optimres$par[j],digits=dig)}
+			optimBootResult[k,"AIC",nGoodBoot] <- round(optimres$AIC,digits=dig)
+			optimBootResult[k,"AICc",nGoodBoot] <- round(optimres$AICc,digits=dig)
+			optimBootResult[k,"BIC",nGoodBoot] <- round(optimres$BIC,digits=dig)
+			optimBootResult[k,"RSS",nGoodBoot] <- round(optimres$value,digits=dig)
+			optimBootResult[k,"R2a",nGoodBoot] <- round(optimres$R2a,digits=dig)
+			optimBootResult[k,"Norm Stat",nGoodBoot] <- round(optimres$normaStat,digits=dig)
+			optimBootResult[k,"Norm p.val",nGoodBoot] <- round(optimres$normaPval,digits=dig)
+			optimBootResult[k,"Pearson",nGoodBoot] <- round(optimres$pearson,digits=dig)
+			optimBootResult[k,"Pea p.val",nGoodBoot] <- round(optimres$pearpval,digits=dig)
+			bootCalculated[k,,nGoodBoot] <- optimres$calculated
+
+			#Fitting validation
+			if (optimBootResult[k,"Norm p.val",nGoodBoot]<0.05 || optimBootResult[k,"Pea p.val",nGoodBoot]<0.05 || length(which(bootCalculated[k,,nGoodBoot]<0)) !=0 ) { flags[k,nGoodBoot]<-"KO"
+			} else {
+				flags[k,nGoodBoot]<-"OK"
+			}#end of if/else on Shap and Corr
+			}#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=T
+	if ( length(which(flags[,nGoodBoot]==0)) != 0 ) badBoot=T
+
+	if (badBoot == F) { 
+			    #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 if badBoot
+}#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)
+nBadBoot = 0
+f=0
+
+for (k in 1:nBoot) {
+
+	cat("Boot Sample : ",k,"\n")
+
+	#Models comparaison
+	#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")}
+		}
+
+	if(verb) cat("Choosen criterion is ",IC,"\n")
+
+	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] <- round(DeltaICvect,digits=dig)
+	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
+
+cat("Bad boot: ",nBadBoot,"\n")
+
+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,choosenModels=choosenModels,IC=IC) 
+
+invisible(res)
+
+} #end of resAverage
+


Property changes on: pkg/R/multiSAR.R
___________________________________________________________________
Name: svn:executable
   + *

Added: pkg/R/rssoptim.R
===================================================================
--- pkg/R/rssoptim.R	                        (rev 0)
+++ pkg/R/rssoptim.R	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,161 @@
+rssoptim <-
+function(model,data,norTest="lillie",graph=T,verb=T,PNGout=NULL){
+
+######################################################
+#                  INPUTS                            #   
+######################################################
+#                                                    #
+# model: the model -list- (ex:  power)               #
+# data: the data -data.frame-                        #
+# norTest : kolmo OR shapiro OR lillie               #
+# graph : print graphics?                            #
+# verb : print info in console?                      #
+######################################################
+
+if (norTest == "lillie") library(nortest)
+
+
+data.name=data$name
+data=data$data[,1:2]
+
+if (verb){
+    cat("**********************************\n")
+    cat("----------------------------------\n")
+    cat("-FITTING: Residual Sum of Squares-\n")
+    cat("< MODEL: ",model$name,">\n")
+    cat("<  DATA: ",data.name,">\n")
+    cat("----------------------------------\n")
+}
+
+    l <- data[[2]]
+    a <- data[[1]]
+
+if (verb){
+    cat("--------------DATAS---------------\n")
+    cat("A:",a,"\n")
+    cat("S:",l,"\n")
+    cat("----------------------------------\n")
+}
+
+	#paramters bounds
+	parLim = model$parLim
+
+	#Transformed initial paramters
+	start <- model$init(data)
+	cat("start :", start,"\n")
+
+
+	for (i in 1:length(start)) { if(parLim[i]!="R"){if(start[i]<=0){ start[i]=0.1 } }  }	
+
+	startMod = transLink(start,parLim)
+    
+	#RSS function
+	rssfun <- model$rssfun
+
+	
+if (verb){
+    cat("------INITIAL VALUES--------------\n")
+    cat(start,"\n")
+    cat("----------------------------------\n")
+    cat("-transformed INITIAL VALUES-------\n")
+    cat(startMod,"\n")
+    cat("----------------------------------\n")
+
+}
+
+
+#if (model$name == "Exponential") res1=optim(start,rssfun,hessian=F,data=data, method="SANN", control=list(maxit=10000)) else 
+res1=optim(startMod,rssfun,hessian=F,data=data,opt=T, method="Nelder-Mead", control=list(maxit=50000)) # CG SANN Nelder-Mead BFGS
+
+#Backtransformation of parameters values
+
+res1$par = backLink(res1$par,parLim)
+
+#if (res1$par[2] <= 0) res1$par[2] <- start[2]
+names(res1$par)=model$paramnames
+
+l.calc=NULL
+l.calc = as.vector(model$fun(res1$par,data))
+residu = as.vector(l - l.calc)
+
+#cat("les residus :" , residu, "\n")
+
+res2 = list(startvalues=start,data=data,model=model,calculated=l.calc,residuals=residu)
+
+#Residuals normality test
+
+normaTest = switch(norTest, "shapiro" = shapiro.test(residu) , "kolmo" = ks.test(residu, "pnorm") , "lillie" = lillie.test(residu) )
+
+#cat("P.Value de " ,norTest, " : " ,normaTest$p,"\n")
+
+#if (norTest == "shapiro") { normaTest <- shapiro.test(residu) } else { normaTest = ks.test(residu, "pnorm") }
+
+#Homogeneity of variance
+cor <- cor.test(residu,data[[1]])
+
+#Nullity of the residuals mean
+
+nullMeanTest <- t.test(residu)
+
+#Calcul des criteres R2a, AIC, AICc, BIC
+
+#variables communes
+n = length(a)
+P = model$paramnumber + 1  # + 1 pour la variance estimee
+
+#R2a
+R2a <- 1 - ((res1$value/n) - P - 1) / (sum((l - mean(l))^2) / (n - 1) )
+
+#AIC
+AIC = n * log(res1$value / n) + 2 * P
+
+#AICc
+AICc = n * log(res1$value / n) + 2*P*(n / (n - P - 1))
+
+#BIC
+BIC = n *log(res1$value / n) + log(n) * P
+
+if(verb){
+        cat("----------FINAL VALUES-----------\n")
+        cat(res1$par,"\n")
+        cat("----------------------------------\n")
+        cat("RSS.value:",res1$value,"\n")
+        cat("----------------------------------\n")
+	cat("------RESIDUALS NORMALITY --------\n")
+	if (norTest == "shapiro") {
+ 	cat("Shapiro Test, W = ",normaTest$statistic,"\n")
+	cat("Shapiro Test, p.value = ",normaTest$p.value,"\n")
+	} else {
+		if (norTest == "kolmo") {
+ 			cat("Kolmogorov Test, D = ",normaTest$statistic,"\n")
+			cat("Kolmogorov Test, p.value = ",normaTest$p.value,"\n")
+		} else {
+			cat("Lilliefors Test, D = ",normaTest$statistic,"\n")
+			cat("Lilliefors Test, p.value = ",normaTest$p.value,"\n")
+		}
+	}
+        cat("------HOMOGENEITY OF VARIANCE ---------\n")
+	cat("Pearson coef. = " ,cor$estimate,"\n")
+	cat("cor. Test p.value  = " ,cor$p.value,"\n")
+	
+	cat("------NULLITY OF RESIDUALS MEAN ---------\n")
+	cat("t = " ,nullMeanTest$statistic,"\n")
+	cat("T Test p.value  = " ,nullMeanTest$p.value,"\n")
+
+
+	cat("--------- CRITERIONS -------------------\n")
+	cat("AIC : ",AIC,"\n")
+	cat("AICc : ",AICc,"\n")
+	cat("BIC : ",BIC,"\n")
+	cat("R2a : ",R2a,"\n")
+	cat("**********************************\n")
+}#end of if verb
+
+    res3 = list(AIC=AIC, AICc=AICc, BIC=BIC, R2a=R2a)
+    res = c(res1,list(pearson=cor$estimate,pearpval=cor$p.value,normaTest=norTest,normaStat=normaTest$statistic,normaPval=normaTest$p.value,tTest=nullMeanTest$statistic,tTestpval=nullMeanTest$p.value),res2,res3,list(data.name=data.name))
+   
+invisible(res)
+
+#END OF RSSOPTIM
+}
+

Added: pkg/R/transLink.R
===================================================================
--- pkg/R/transLink.R	                        (rev 0)
+++ pkg/R/transLink.R	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,30 @@
+transLink <-
+function(x,boundType) {
+
+	# fonction logit
+	logit=function(y) {
+		log(y/(1-y))
+	}#end of logit
+
+if (length(x) ==1) {
+	res = switch(boundType,
+		R = x,
+		Rplus = log(x),
+		unif = logit(x)
+	)#end of switch
+} else {
+
+res=vector()
+
+	for (i in 1:length(x)) {
+
+		res[i] = switch(boundType[i],
+				R = x[i],
+				Rplus = log(x[i]),
+				unif = logit(x[i])
+			)#end of switch
+	}#end of for
+res
+}#end of if/else
+}#end of transLink
+

Added: pkg/data/asymp.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/asymp.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/data.arr.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/data.arr.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/data.atl.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/data.atl.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/data.galap.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/data.galap.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/data.glea.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/data.glea.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/data.gleas.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/data.gleas.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/expo.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/expo.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/logist.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/logist.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/lomolino.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/lomolino.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/monod.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/monod.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/negexpo.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/negexpo.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/power.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/power.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/ratio.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/ratio.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/data/weibull.rda
===================================================================
(Binary files differ)


Property changes on: pkg/data/weibull.rda
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: pkg/man/asymp.Rd
===================================================================
--- pkg/man/asymp.Rd	                        (rev 0)
+++ pkg/man/asymp.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,41 @@
+\name{asymp}
+\Rdversion{1.1}
+\alias{asymp}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(asymp)}
+\format{
+  The format is:
+List of 9
+ $ name       : chr "asymptotic regression"
+ $ formula    :  expression(S == c - z * f^A)
+ $ paramnumber: num 3
+ $ paramnames : chr [1:3] "c" "z" "f"
+ $ limits     : num [1:2, 1:3] 0 Inf 0 Inf 0 ...
+ $ parLim     : chr [1:3] "Rplus" "Rplus" "unif"
+ $ fun        :function (par, data)  
+  ..- attr(*, "source")= chr "function(par,data){if(length(data)>1) d=data[[1]] else d=data; s = par[1] - par[2] * par[3]^data[[1]]; as.vector(s)}"
+ $ rssfun     :function (par, data, opt)  
+  ..- attr(*, "source")= chr "function(par,data,opt){if(opt)par=backLink(par,asymp$parLim) ; sum( ( data[[2]] - (par[1] - par[2] * par[3]^data[[1]]) )^2  ) }"| __truncated__
+ $ init       :function (data)  
+  ..- attr(*, "source")= chr [1:10] "function(data){#Ratkowsky 1983 p178" ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(asymp)
+## maybe str(asymp) ; plot(asymp) ...
+}
+\keyword{datasets}

Added: pkg/man/data.arr.Rd
===================================================================
--- pkg/man/data.arr.Rd	                        (rev 0)
+++ pkg/man/data.arr.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,33 @@
+\name{data.arr}
+\Rdversion{1.1}
+\alias{data.arr}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(data.arr)}
+\format{
+  The format is:
+List of 2
+ $ name: chr "Arrhenius 1921"
+ $ data:'data.frame':	10 obs. of  2 variables:
+  ..$ a: num [1:10] 1 2 4 8 16 32 64 128 256 300
+  ..$ s: num [1:10] 1.4 2 2.6 3.7 5.2 6.8 8.8 10.3 14.5 16
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(data.arr)
+## maybe str(data.arr) ; plot(data.arr) ...
+}
+\keyword{datasets}

Added: pkg/man/data.atl.Rd
===================================================================
--- pkg/man/data.atl.Rd	                        (rev 0)
+++ pkg/man/data.atl.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,33 @@
+\name{data.atl}
+\Rdversion{1.1}
+\alias{data.atl}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(data.atl)}
+\format{
+  The format is:
+List of 2
+ $ name: chr "Preston 1962 (atlantic)"
+ $ data:'data.frame':	11 obs. of  2 variables:
+  ..$ a: num [1:11] 120 150 233 304 465 ...
+  ..$ s: num [1:11] 29 35 35 36 35 37 79 99 74 124 ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(data.atl)
+## maybe str(data.atl) ; plot(data.atl) ...
+}
+\keyword{datasets}

Added: pkg/man/data.galap.Rd
===================================================================
--- pkg/man/data.galap.Rd	                        (rev 0)
+++ pkg/man/data.galap.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,33 @@
+\name{data.galap}
+\Rdversion{1.1}
+\alias{data.galap}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(data.galap)}
+\format{
+  The format is:
+List of 2
+ $ name: chr "Preston 1962 (galapagos)"
+ $ data:'data.frame':	16 obs. of  2 variables:
+  ..$ a: num [1:16] 0.2 0.9 1 1.8 1.87 4.4 7.1 7.5 18 20 ...
+  ..$ s: num [1:16] 48 7 52 14 42 22 103 48 79 119 ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(data.galap)
+## maybe str(data.galap) ; plot(data.galap) ...
+}
+\keyword{datasets}

Added: pkg/man/data.glea.Rd
===================================================================
--- pkg/man/data.glea.Rd	                        (rev 0)
+++ pkg/man/data.glea.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,33 @@
+\name{data.glea}
+\Rdversion{1.1}
+\alias{data.glea}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(data.glea)}
+\format{
+  The format is:
+List of 2
+ $ name: chr "gleason"
+ $ data:'data.frame':	19 obs. of  2 variables:
+  ..$ a: num [1:19] 1 2 3 4 5 6 8 10 12 15 ...
+  ..$ s: num [1:19] 4.38 5.82 6.9 7.6 8.21 ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(data.glea)
+## maybe str(data.glea) ; plot(data.glea) ...
+}
+\keyword{datasets}

Added: pkg/man/data.gleas.Rd
===================================================================
--- pkg/man/data.gleas.Rd	                        (rev 0)
+++ pkg/man/data.gleas.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,33 @@
+\name{data.gleas}
+\Rdversion{1.1}
+\alias{data.gleas}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(data.gleas)}
+\format{
+  The format is:
+List of 2
+ $ name: chr "gleason"
+ $ data:'data.frame':	19 obs. of  2 variables:
+  ..$ a: num [1:19] 1 2 3 4 5 6 8 10 12 15 ...
+  ..$ s: num [1:19] 4.38 6.67 8.75 10.37 11.29 ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(data.gleas)
+## maybe str(data.gleas) ; plot(data.gleas) ...
+}
+\keyword{datasets}

Added: pkg/man/expo.Rd
===================================================================
--- pkg/man/expo.Rd	                        (rev 0)
+++ pkg/man/expo.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,41 @@
+\name{expo}
+\Rdversion{1.1}
+\alias{expo}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(expo)}
+\format{
+  The format is:
+List of 9
+ $ name       : chr "Exponential"
+ $ formula    :  expression(s == z * log(a) + c)
+ $ paramnumber: num 2
+ $ paramnames : chr [1:2] "c" "z"
+ $ limits     : num [1:2, 1:2] 0 Inf 0 Inf
+ $ parLim     : chr [1:2] "R" "Rplus"
+ $ fun        :function (par, data)  
+  ..- attr(*, "source")= chr "function(par,data){if(length(data)>1) d=data[[1]] else d=data; s = par[2] * log(d) + par[1]; names(s)=c(\"s.expo\"); as.vector("| __truncated__
+ $ rssfun     :function (par, data, opt)  
+  ..- attr(*, "source")= chr "function(par,data,opt){if(opt)par=backLink(par,expo$parLim) ; sum(   (data[[2]] - (par[2] * log(data[[1]]) + par[1]) ) ^2     )"| __truncated__
+ $ init       :function (data)  
+  ..- attr(*, "source")= chr [1:7] "function(data){" ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(expo)
+## maybe str(expo) ; plot(expo) ...
+}
+\keyword{datasets}

Added: pkg/man/logist.Rd
===================================================================
--- pkg/man/logist.Rd	                        (rev 0)
+++ pkg/man/logist.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,43 @@
+\name{logist}
+\Rdversion{1.1}
+\alias{logist}
+\docType{data}
+\title{
+%%   ~~ data name/kind ... ~~
+}
+\description{
+%%  ~~ A concise (1-5 lines) description of the dataset. ~~
+}
+\usage{data(logist)}
+\format{
+  The format is:
+List of 10
+ $ name       : chr "Logistic function"
+ $ formula    :  expression(S == over(c, (1 + exp(-z * A + f))))
+ $ Formula    :Class 'formula' length 3 s ~ over(par1, (1 + exp(-par2 * a + par3)))
+  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
+ $ paramnumber: num 3
+ $ paramnames : chr [1:3] "c" "z" "f"
+ $ limits     : num [1:2, 1:3] 0 Inf 0 Inf -Inf ...
+ $ parLim     : chr [1:3] "Rplus" "Rplus" "R"
+ $ fun        :function (par, data)  
+  ..- attr(*, "source")= chr "function(par,data){if(length(data)>1) d=data[[1]] else d=data; s = par[1] / (1 + exp(-par[2]*d+par[3])) ; as.vector(s)}"
+ $ rssfun     :function (par, data, opt)  
+  ..- attr(*, "source")= chr "function(par,data,opt){if(opt)par=backLink(par,logist$parLim) ; sum( ( data[[2]] - (par[1] / (1 + exp(-par[2]*data[[1]]+par[3])"| __truncated__
+ $ init       :function (data)  
+  ..- attr(*, "source")= chr [1:14] "function(data){" ...
+}
+\details{
+%%  ~~ If necessary, more details than the __description__ above ~~
+}
+\source{
+%%  ~~ reference to a publication or URL from which the data were obtained ~~
+}
+\references{
+%%  ~~ possibly secondary sources and usages ~~
+}
+\examples{
+data(logist)
+## maybe str(logist) ; plot(logist) ...
+}
+\keyword{datasets}

Added: pkg/man/lomolino.Rd
===================================================================
--- pkg/man/lomolino.Rd	                        (rev 0)
+++ pkg/man/lomolino.Rd	2009-09-22 12:32:47 UTC (rev 9)
@@ -0,0 +1,41 @@
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/mmsar -r 9


More information about the Mmsar-commits mailing list