[Mmsar-commits] r3 - in pkg: . pkg pkg/R pkg/data pkg/man www
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 15 10:12:37 CEST 2009
Author: fguilhaumon
Date: 2009-09-15 10:12:35 +0200 (Tue, 15 Sep 2009)
New Revision: 3
Added:
pkg/README
pkg/pkg/
pkg/pkg/DESCRIPTION
pkg/pkg/R/
pkg/pkg/R/backLink.R
pkg/pkg/R/multiSAR.R
pkg/pkg/R/rssoptim.R
pkg/pkg/R/transLink.R
pkg/pkg/Read-and-delete-me
pkg/pkg/data/
pkg/pkg/data/asymp.rda
pkg/pkg/data/data.arr.rda
pkg/pkg/data/data.atl.rda
pkg/pkg/data/data.galap.rda
pkg/pkg/data/data.glea.rda
pkg/pkg/data/data.gleas.rda
pkg/pkg/data/expo.rda
pkg/pkg/data/logist.rda
pkg/pkg/data/lomolino.rda
pkg/pkg/data/monod.rda
pkg/pkg/data/negexpo.rda
pkg/pkg/data/power.rda
pkg/pkg/data/ratio.rda
pkg/pkg/data/weibull.rda
pkg/pkg/man/
pkg/pkg/man/asymp.Rd
pkg/pkg/man/data.arr.Rd
pkg/pkg/man/data.atl.Rd
pkg/pkg/man/data.galap.Rd
pkg/pkg/man/data.glea.Rd
pkg/pkg/man/data.gleas.Rd
pkg/pkg/man/expo.Rd
pkg/pkg/man/logist.Rd
pkg/pkg/man/lomolino.Rd
pkg/pkg/man/mmSAR.Rd
pkg/pkg/man/monod.Rd
pkg/pkg/man/multiSAR.Rd
pkg/pkg/man/negexpo.Rd
pkg/pkg/man/power.Rd
pkg/pkg/man/ratio.Rd
pkg/pkg/man/rssoptim.Rd
pkg/pkg/man/weibull.Rd
pkg/www/
pkg/www/index.php
Log:
Added: pkg/README
===================================================================
--- pkg/README (rev 0)
+++ pkg/README 2009-09-15 08:12:35 UTC (rev 3)
@@ -0,0 +1,58 @@
+ R-Forge SVN README
+
+
+(See "http://download.r-forge.r-project.org/manuals/R-Forge_Manual.pdf"
+ for detailed information on registering a new project.)
+
+1. Introduction
+-----------------------------------------------------------------------
+R is free software distributed under a GNU-style copyleft. R-Forge is
+a central platform for the development of R packages, R-related
+software and further projects. Among many other web-based features it
+provides facilities for collaborative source code management via
+Subversion (SVN).
+
+2. The directory you're in
+-----------------------------------------------------------------------
+This is the repository of your project. It contains two important
+pre-defined directories namely 'pkg' and 'www'. These directories must
+not be deleted otherwise R-Forge's core functionality will not be
+available (i.e., daily checking and building of your package or the
+project websites).
+'pkg' and 'www' are standardized and therefore are going to be
+described in this README. The rest of your repository can be used as
+you like.
+
+3. 'pkg' directory
+-----------------------------------------------------------------------
+To make use of the package building and checking feature the package
+source code has to be put into the 'pkg' directory of your repository
+(i.e., 'pkg/DESCRIPTION', 'pkg/R', 'pkg/man', etc.) or, alternatively,
+a subdirectory of 'pkg'. The latter structure allows for having more
+than one package in a single project, e.g., if a project consists of
+the packages foo and bar then the source code will be located in
+'pkg/foo' and 'pkg/bar', respectively.
+
+R-Forge automatically examines the 'pkg' directory of every repository
+and builds the package sources as well as the package binaries on a
+daily basis for Mac OSX and Windows (if applicable). The package builds
+are provided in the 'R Packages' tab for download or can be installed
+directly in R from a CRAN-style repository using
+'install.packages("foo", repos="http://R-Forge.R-project.org")'.
+Furthermore, in the 'R Packages' tab developers can examine logs
+generated by the build and check process on different platforms.
+
+4. 'www' directory
+-----------------------------------------------------------------------
+Developers may present their work on a subdomain of R-Forge, e.g.,
+'http://foo.R-Forge.R-project.org', or via a link to an external
+website.
+
+This directory contains the project homepage which gets updated hourly
+on R-Forge, so please take into consideration that it will not be
+available right after you commit your changes or additions.
+
+5. Help
+-----------------------------------------------------------------------
+If you need help don't hesitate to contact us at R-Forge at R-project.org
+
Added: pkg/pkg/DESCRIPTION
===================================================================
--- pkg/pkg/DESCRIPTION (rev 0)
+++ pkg/pkg/DESCRIPTION 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/R/backLink.R
===================================================================
--- pkg/pkg/R/backLink.R (rev 0)
+++ pkg/pkg/R/backLink.R 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/R/multiSAR.R
===================================================================
--- pkg/pkg/R/multiSAR.R (rev 0)
+++ pkg/pkg/R/multiSAR.R 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/R/multiSAR.R
___________________________________________________________________
Name: svn:executable
+
Added: pkg/pkg/R/rssoptim.R
===================================================================
--- pkg/pkg/R/rssoptim.R (rev 0)
+++ pkg/pkg/R/rssoptim.R 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/R/transLink.R
===================================================================
--- pkg/pkg/R/transLink.R (rev 0)
+++ pkg/pkg/R/transLink.R 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/Read-and-delete-me
===================================================================
--- pkg/pkg/Read-and-delete-me (rev 0)
+++ pkg/pkg/Read-and-delete-me 2009-09-15 08:12:35 UTC (rev 3)
@@ -0,0 +1,9 @@
+* Edit the help file skeletons in 'man', possibly combining help files
+ for multiple functions.
+* Put any C/C++/Fortran code in 'src'.
+* If you have compiled code, add a .First.lib() function in 'R' to load
+ the shared library.
+* Run R CMD build to build the package tarball.
+* Run R CMD check to check the package tarball.
+
+Read "Writing R Extensions" for more information.
Added: pkg/pkg/data/asymp.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/asymp.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/data.arr.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/data.arr.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/data.atl.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/data.atl.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/data.galap.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/data.galap.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/data.glea.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/data.glea.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/data.gleas.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/data.gleas.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/expo.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/expo.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/logist.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/logist.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/lomolino.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/lomolino.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/monod.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/monod.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/negexpo.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/negexpo.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/power.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/power.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/ratio.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/ratio.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/data/weibull.rda
===================================================================
(Binary files differ)
Property changes on: pkg/pkg/data/weibull.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/pkg/man/asymp.Rd
===================================================================
--- pkg/pkg/man/asymp.Rd (rev 0)
+++ pkg/pkg/man/asymp.Rd 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/man/data.arr.Rd
===================================================================
--- pkg/pkg/man/data.arr.Rd (rev 0)
+++ pkg/pkg/man/data.arr.Rd 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/man/data.atl.Rd
===================================================================
--- pkg/pkg/man/data.atl.Rd (rev 0)
+++ pkg/pkg/man/data.atl.Rd 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/man/data.galap.Rd
===================================================================
--- pkg/pkg/man/data.galap.Rd (rev 0)
+++ pkg/pkg/man/data.galap.Rd 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/man/data.glea.Rd
===================================================================
--- pkg/pkg/man/data.glea.Rd (rev 0)
+++ pkg/pkg/man/data.glea.Rd 2009-09-15 08:12:35 UTC (rev 3)
@@ -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/pkg/man/data.gleas.Rd
===================================================================
--- pkg/pkg/man/data.gleas.Rd (rev 0)
+++ pkg/pkg/man/data.gleas.Rd 2009-09-15 08:12:35 UTC (rev 3)
@@ -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:
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/mmsar -r 3
More information about the Mmsar-commits
mailing list