[Mmsar-commits] r13 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 6 22:24:25 CET 2010
Author: fguilhaumon
Date: 2010-02-06 22:24:25 +0100 (Sat, 06 Feb 2010)
New Revision: 13
Modified:
pkg/DESCRIPTION
pkg/R/multiSAR.R
pkg/R/rssoptim.R
pkg/man/rssoptim.Rd
Log:
update manpage for rssoptim
code clean for other functions
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-02-06 19:25:04 UTC (rev 12)
+++ pkg/DESCRIPTION 2010-02-06 21:24:25 UTC (rev 13)
@@ -5,6 +5,7 @@
Date: 2009-09-01
Author: François Guilhaumon
Maintainer: François Guilhaumon <francois.guilhaumon at univ-montp2.fr>
+Depends: R (>= 2.5.0), numDeriv, nortest
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
Modified: pkg/R/multiSAR.R
===================================================================
--- pkg/R/multiSAR.R 2010-02-06 19:25:04 UTC (rev 12)
+++ pkg/R/multiSAR.R 2010-02-06 21:24:25 UTC (rev 13)
@@ -1,6 +1,9 @@
multiSAR <-
-function(modelList=c("power","expo","negexpo","monod","logist","ratio","lomolino","weibull"),data=gleason,nBoot=999,verb=T,crit="Bayes") {
+function(modelList=c("power","expo","negexpo","monod","logist","ratio","lomolino","weibull"),data,nBoot=999,verb=F,crit="Bayes") {
+library(numDeriv)
+
+
############Criteria must be "Info" for AIC/AICc or "Bayes" for BIC
nlig <- length(modelList)
Modified: pkg/R/rssoptim.R
===================================================================
--- pkg/R/rssoptim.R 2010-02-06 19:25:04 UTC (rev 12)
+++ pkg/R/rssoptim.R 2010-02-06 21:24:25 UTC (rev 13)
@@ -1,5 +1,5 @@
rssoptim <-
-function(model,data,norTest="lillie",graph=T,verb=T,PNGout=NULL){
+function(model,data,norTest="lillie",verb=T){
######################################################
# INPUTS #
@@ -7,35 +7,34 @@
# #
# model: the model -list- (ex: power) #
# data: the data -data.frame- #
-# norTest : kolmo OR shapiro OR lillie #
-# graph : print graphics? #
+# norTest : shapiro OR lillie #
# verb : print info in console? #
######################################################
-if (norTest == "lillie") library(nortest)
+ if (norTest == "lillie") library(nortest)
-data.name=data$name
-data=data$data[,1:2]
+ 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")
-}
+ 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]]
+ l <- data[[2]]
+ a <- data[[1]]
-if (verb){
- cat("--------------DATAS---------------\n")
- cat("A:",a,"\n")
- cat("S:",l,"\n")
- cat("----------------------------------\n")
-}
+ if (verb){
+ cat("--------------DATAS---------------\n")
+ cat("A:",a,"\n")
+ cat("S:",l,"\n")
+ cat("----------------------------------\n")
+ }
#paramters bounds
parLim = model$parLim
@@ -44,7 +43,6 @@
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)
@@ -53,106 +51,86 @@
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 (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
+ 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
+ #Backtransformation of parameters values
-res1$par = backLink(res1$par,parLim)
+ res1$par = backLink(res1$par,parLim)
-#if (res1$par[2] <= 0) res1$par[2] <- start[2]
-names(res1$par)=model$paramnames
+ #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)
+ 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)
-res2 = list(startvalues=start,data=data,model=model,calculated=l.calc,residuals=residu)
+ #Residuals normality test
-#Residuals normality test
+ normaTest = switch(norTest, "shapiro" = shapiro.test(residu) , "lillie" = lillie.test(residu) )
-normaTest = switch(norTest, "shapiro" = shapiro.test(residu) , "kolmo" = ks.test(residu, "pnorm") , "lillie" = lillie.test(residu) )
+ #Homogeneity of variance
+ cor <- cor.test(residu,data[[1]])
-#cat("P.Value de " ,norTest, " : " ,normaTest$p,"\n")
+ #Calcul des criteres R2a, AIC, AICc, BIC
-#if (norTest == "shapiro") { normaTest <- shapiro.test(residu) } else { normaTest = ks.test(residu, "pnorm") }
+ #variables communes
+ n = length(a)
+ P = model$paramnumber + 1 # + 1 pour la variance estimee
-#Homogeneity of variance
-cor <- cor.test(residu,data[[1]])
+ #R2 (this definition of the R2 alloaws the comparison of all models)
+ R2 <- 1 - ( (res1$value) / sum((l - mean(l))^2) )
-#Nullity of the residuals mean
+ #AIC
+ AIC = n * log(res1$value / n) + 2 * P
-nullMeanTest <- t.test(residu)
+ #AICc
+ AICc = n * log(res1$value / n) + 2*P*(n / (n - P - 1))
-#Calcul des criteres R2a, AIC, AICc, BIC
+ #BIC
+ BIC = n *log(res1$value / n) + log(n) * P
-#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")
+ 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 {
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("------HOMOGENEITY OF VARIANCE ---------\n")
+ cat("Pearson coef. = " ,cor$estimate,"\n")
+ cat("cor. Test p.value = " ,cor$p.value,"\n")
+
+ cat("--------- CRITERIA -------------------\n")
+ cat("AIC : ",AIC,"\n")
+ cat("AICc : ",AICc,"\n")
+ cat("BIC : ",BIC,"\n")
+ cat("R2 : ",R2,"\n")
+ cat("**********************************\n")
+ }#end of if verb
-
- 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))
+ res3 = list(AIC=AIC, AICc=AICc, BIC=BIC, R2=R2)
+ res = c(res1,list(pearson=cor$estimate,pearpval=cor$p.value,normaTest=norTest,normaStat=normaTest$statistic,normaPval=normaTest$p.value,res2,res3,list(data.name=data.name))
invisible(res)
Modified: pkg/man/rssoptim.Rd
===================================================================
--- pkg/man/rssoptim.Rd 2010-02-06 19:25:04 UTC (rev 12)
+++ pkg/man/rssoptim.Rd 2010-02-06 21:24:25 UTC (rev 13)
@@ -1,229 +1,90 @@
\name{rssoptim}
\Rdversion{1.1}
\alias{rssoptim}
-%- Also NEED an '\alias' for EACH other topic documented here.
-\title{
-%% ~~function to do ... ~~
+\title{The Residual Sum of Squares optimisation function}
+\description{The rssoptim function is used to fit SAR models (as mmSAR model objects) to dataset (as mmSAR data objects). This function allow the optimisation of the model parameters, fitted values calculation, information criteria (AIC, AICc, BIC) and R² calculations and fit evaluation with residuals normality and homoscedasticity tests.
}
-\description{
-%% ~~ A concise (1-5 lines) description of what the function does. ~~
-}
\usage{
-rssoptim(model, data, norTest = "lillie", graph = T, verb = T, PNGout = NULL)
+rssoptim(model, data, norTest = "lillie", verb = T)
}
-%- maybe also 'usage' for other objects documented here.
\arguments{
- \item{model}{
-%% ~~Describe \code{model} here~~
+ \item{model}{A mmSAR model object representing the model to fit.
}
- \item{data}{
-%% ~~Describe \code{data} here~~
+ \item{data}{A mmSAR data object representing the data to fit the model on.
}
- \item{norTest}{
-%% ~~Describe \code{norTest} here~~
+ \item{norTest}{A character string from "lillie" or "shapiro" for respectively test for the normality of the residuals with a Lilliefors or Shapiro test.
}
- \item{graph}{
-%% ~~Describe \code{graph} here~~
+ \item{verb}{Bolean : should the function reports informations while running ?
}
- \item{verb}{
-%% ~~Describe \code{verb} here~~
}
- \item{PNGout}{
-%% ~~Describe \code{PNGout} here~~
-}
-}
\details{
-%% ~~ If necessary, more details than the description above ~~
+The rssoptim function is basically a wrapper for the optim R built-in function. rssoptim obtains model parameters estimates by minimizing the residual sum of squares using the unconstrained Nelder–Mead optimization algorithm. Assuming normality of the observations, this approach produces optimal maximum likelihood estimates of model parameters. rssoptim provides two tests for the normality of the residuals : the Lilliefors extension of the Kolmogorov normality test, which is advocated when sample size is large or when the data show a substantial variability (e.g. continental scale studies) and the Shapiro-Wilk test for normality, which focuses on skewness and kurtosis of the empirical distribution of the residuals and is useful for small sample size or when data results from small scale sampling. mmSAR tests for homoscedasticity by evaluating the correlation between residuals magnitude and areas (Pearson’s product moment correlation coefficient).
}
\value{
-%% ~Describe the value returned
-%% If it is a LIST, use
-%% \item{comp1 }{Description of 'comp1'}
-%% \item{comp2 }{Description of 'comp2'}
-%% ...
-}
-\references{
-%% ~put references to the literature/web site here ~
-}
-\author{
-%% ~~who you are~~
-}
-\note{
-%% ~~further notes~~
-}
+A list :
-%% ~Make other sections like Warning with \section{Warning }{....} ~
+$par : the parameters estimates.
-\seealso{
-%% ~~objects to See Also as \code{\link{help}}, ~~~
-}
-\examples{
-##---- Should be DIRECTLY executable !! ----
-##-- ==> Define data, use random,
-##-- or do help(data=index) for the standard data sets.
+$value : the value of the Residual Sum of Squares (RSS).
-## The function is currently defined as
-function(model,data,norTest="lillie",graph=T,verb=T,PNGout=NULL){
+$counts : number of iterations for the convergence of the fitting algorithm
-######################################################
-# 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? #
-######################################################
+$convergence : optim convergence code (only 0 is OK).
-if (norTest == "lillie") library(nortest)
+$message : optim convergence message (only NULL is OK).
+$pearson : Pearson’s product moment correlation coefficient value.
-data.name=data$name
-data=data$data[,1:2]
+$pearpval : Pearson’s product moment correlation coefficient p.value.
-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")
- }
+$normaTest : residuals normality test used.
- l <- data[[2]]
- a <- data[[1]]
+$normaStat : residuals normality test statistic value.
-if (verb){
- cat("--------------DATAS---------------\n")
- cat("A:",a,"\n")
- cat("S:",l,"\n")
- cat("----------------------------------\n")
- }
+$normaPval : residuals normality test statistic p.value.
- #paramters bounds
- parLim = model$parLim
+$startvalues : starting values for the fitting algorithm.
- #Transformed initial paramters
- start <- model$init(data)
- cat("start :", start,"\n")
+$data : the data.
+$model : the mmSAR model object.
- for (i in 1:length(start)) { if(parLim[i]!="R"){if(start[i]<=0){ start[i]=0.1 } } }
+$calculated : vector of fitted values.
- startMod = transLink(start,parLim)
-
- #RSS function
- rssfun <- model$rssfun
+$residuals : vector of residuals.
-
-if (verb){
- cat("------INITIAL VALUES--------------\n")
- cat(start,"\n")
- cat("----------------------------------\n")
- cat("-transformed INITIAL VALUES-------\n")
- cat(startMod,"\n")
- cat("----------------------------------\n")
+$AIC : AIC value.
- }
+$AICc : AICc value.
+$BIC : BIC value.
-#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
+$R2a : R² value.
-#Backtransformation of parameters values
+$data.name : the name of the dataset (from the mmSAR data object).
-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)
+\author{francois Guilhaumon
+}
-#cat("les residus :" , residu, "\n")
+\seealso{
+}
-res2 = list(startvalues=start,data=data,model=model,calculated=l.calc,residuals=residu)
+\examples{
-#Residuals normality test
+#FITTING THE POWER MODEL TO THE GALAPAGOS DATASET
-normaTest = switch(norTest, "shapiro" = shapiro.test(residu) , "kolmo" = ks.test(residu, "pnorm") , "lillie" = lillie.test(residu) )
+#loading the data
+data(data.galap)
-#cat("P.Value de " ,norTest, " : " ,normaTest$p,"\n")
+#loading the model
+data(power)
-#if (norTest == "shapiro") { normaTest <- shapiro.test(residu) } else { normaTest = ks.test(residu, "pnorm") }
+#fitting the model with verbosity
+rssoptim(power,data.galap,verb=T)
-#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
- }
}
-% Add one or more standard keywords, see file 'KEYWORDS' in the
-% R documentation directory.
-\keyword{ ~kwd1 }
-\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
More information about the Mmsar-commits
mailing list