Merge pull request #85 from paternogbc/version_0.3.0

Version 0.3.0
paternogbc · Nov 16, 2015 · 168705e · 168705e
2 parents 0ba866a + 35aee66
commit 168705e
Show file tree

Hide file tree

Showing 28 changed files with 162 additions and 544 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -1,8 +1,9 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 NEWS\.md
-^developing$
+devloping
 
 ^\.travis\.yml$
 ^data-raw$
 CODE_OF_CONDUCT.md
+
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: sensiPhy
 Title: Sensitivity Analysis for Comparative Methods
-Version: 0.2.2
-Date: 2015-06-08
+Version: 0.3.0
+Date: 2015-9-11
 Authors@R: c(
     person("Gustavo","Paterno", email = "[email protected]", role = c("cre","aut")),
     person("Gijsbert", "Werner", email = "[email protected]", role = "aut"),
@@ -21,4 +21,4 @@ Imports:
 Suggests:
     knitr (>= 1.10),
     testthat
-VignetteBuilder: knitr
+VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
@@ -14,8 +14,8 @@ export(clade_phyglm)
 export(clade_phylm)
 export(influ_phyglm)
 export(influ_phylm)
+export(intra_phyglm)
 export(intra_phylm)
-export(intra_phyloglm)
 export(match_dataphy)
 export(samp_phyglm)
 export(samp_phylm)

diff --git a/R/clade_phylolm.R b/R/clade_phylolm.R
@@ -52,25 +52,6 @@
 #' reported.
 #' @return \code{data}: Original full dataset.
 #' @return \code{errors}: Clades where deletion resulted in errors.
-#' @examples
-#' \dontrun{
-#' library(sensiPhy);library(phylolm)
-#'
-#' #Generate a random tree
-#' set.seed(2468)
-#' tree <- rtree(100)
-#'
-#' #Generate random predictor variable (pred), evolving according to a BM model.
-#' pred<- rTraitCont(tree,root.value=0,sigma=1,model="BM")
-#'
-#' #Generate one continous traits
-#' cont_trait <- pred + rTraitCont(tree,model="BM",sigma=2.5)
-#' fam <- rep(c("fam1","fam2","fam3","fam4","fam5"),each=20)
-#' dat<-data.frame(pred,cont_trait,fam)
-#'
-#' #Determine influential clades:
-#' clade.test <- clade_phylm(cont_trait~pred,data=dat,phy=tree,clade.col="fam")
-#' }
 #' @author Gustavo Paterno
 #' @seealso \code{\link[phylolm]{phylolm}}, \code{\link{influ_phylm}},
 #' \code{\link{sensi_plot}}

diff --git a/R/influ_phyloglm.R b/R/influ_phyloglm.R
@@ -54,37 +54,6 @@
 #' (i.e. \code{alpha}) are reported.
 #' @return \code{data}: Original full dataset.
 #' @return \code{errors}: Species where deletion resulted in errors.
-#' @examples
-#' \dontrun{
-#' library(sensiPhy)
-#'
-#' #Generate a random tree
-#' set.seed(2468)
-#' tree <- rtree(100)
-#'
-#' #Generate random predictor variable (pred), evolving according to a BM model.
-#' pred<- rTraitCont(tree,root.value=0,sigma=1,model="BM")
-#'
-#' #Generate two continous traits, one evolving highly correlated with the
-#' #predictor (trait 1), and one evolving more randomly (trait 2)
-#' cont_trait1 <- pred + rTraitCont(tree,model="BM",sigma=0.1)
-#' cont_trait2 <- pred + rTraitCont(tree,model="BM",sigma=10)
-#'
-#' #Generate two binary traits, one highly correlated to pred (trait 1), the other less.
-#' bin_trait1<-rbinTrait(n=1,tree,beta=c(-1,0.5),alpha=0.1,
-#'                      X=cbind(rep(1,length(tree$tip.label)),pred))
-#' bin_trait2<-rbinTrait(n=1,tree,beta=c(-1,0.5),alpha=5,
-#'                       X=cbind(rep(1,length(tree$tip.label)),pred))
-#' dat<-data.frame(pred,cont_trait1,cont_trait2,bin_trait1,bin_trait2)
-#'
-#' #Determine influential species for both regressions.
-#' fit1<-influ_phyglm(bin_trait1~pred,data = dat,phy = tree)
-#' fit2<-influ_phyglm(bin_trait2~pred,data = dat,phy = tree)
-#'
-#' #For purposes of comparison the full model output from phylolm:
-#' summary(phyloglm(bin_trait1~pred,data = dat,phy = tree,method = "logistic_MPLE"),btol=50)
-#' summary(phyloglm(bin_trait2~pred,data = dat,phy = tree,method = "logistic_MPLE"),btol=50)
-#' }
 #' @author Gustavo Paterno & Gijsbert D.A. Werner
 #' @seealso \code{\link[phylolm]{phyloglm}}, \code{\link{samp_phyglm}},
 #' \code{\link{influ_phylm}}, \code{\link{sensi_plot}}
@@ -150,7 +119,8 @@ influ_phyglm <- function(formula,data,phy,btol=50,cutoff=2,track=TRUE,...){
                         aic.mod              <- mod$aic
                         optpar               <- mod$alpha
 
-                        if(track==TRUE) (print(paste(i," / ",N,sep="")))
+                        if(track==TRUE)
+                            cat("\r", i," / ", N)
 
                         #Stores values for eacht simulation
                         estim.simu <- data.frame(sp, intercept, DFintercept, intercept.perc,

diff --git a/R/influ_phylolm.R b/R/influ_phylolm.R
@@ -56,24 +56,6 @@
 #' reported.
 #' @return \code{data}: Original full dataset.
 #' @return \code{errors}: Species where deletion resulted in errors.
-#' @examples
-#' 
-#' library(sensiPhy)
-#' 
-#' # Loading data and phylogeny:
-#' data(alien)
-#' trait <- alien$data
-#' phy <- alien$phy[[1]]
-#' 
-#' # Run sensitivity analysis (influential species)
-#' out <- influ_phylm( log10(mass) ~ log10(georange), data = trait, phy = phy) 
-#'                      
-#' # Check most influential species for slope and intercept:
-#' summary(out)
-#' 
-#' # Plot results:
-#' sensi_plot(out)
-
 #' @author Gustavo Paterno & Gijsbert D.A. Werner
 #' @seealso \code{\link[phylolm]{phylolm}}, \code{\link{samp_phylm}},
 #' \code{\link{influ_phylm}},\code{\link{sensi_plot}}
@@ -141,7 +123,8 @@ influ_phylm <- function(formula,data,phy,model="lambda",cutoff=2,track=TRUE,...)
                             optpar               <- mod$optpar
                         }
 
-                        if(track==TRUE) (print(paste(i," / ",N,sep="")))
+                        if(track==TRUE)
+                        cat("\r", i," / ", N)
 
                         # Stores values for each simulation
                         estim.simu <- data.frame(sp, intercept, DFintercept, intercept.perc,

diff --git a/R/intra_phyloglm.R b/R/intra_phyloglm.R
@@ -8,6 +8,8 @@
 #' @param phy A phylogeny (class 'phylo', see ?\code{ape}).
 #' @param Vx Name of the column containing the standard deviation or the standard error of the predictor 
 #' variable. When information is not available for one taxon, the value can be 0 or \code{NA}
+#' @param Vy Name of the column containing the standard deviation or the standard error of the response 
+#' variable. When information is not available for one taxon, the value can be 0 or \code{NA}
 #' @param times Number of times to repeat the analysis generating a random value for the predictor variable.
 #' If NULL, \code{times} = 2
 #' @param distrib A character string indicating which distribution to use to generate a random value for the response 
@@ -43,28 +45,21 @@
 #' @return \code{stats}: Statistics for model parameters. \code{sd_intra} is the standard deviation 
 #' due to intraspecific variation. \code{CI_low} and \code{CI_high} are the lower and upper limits 
 #' of the 95% confidence interval.
-#' @examples
-#' \dontrun{
-#' library(sensiPhy)
-#' }
-#' 
 #' @author Caterina Penone & Pablo Ariel Martinez
 #' @seealso \code{\link{sensi_plot}}
 #' @references Here still: reference to phylolm paper + our own?
 #' @export
 
 
-intra_phyloglm <- function(formula,data,phy,
-                          Vx=NULL,times=2,
-                          distrib="uniform",btol=50,track=TRUE,...){
+intra_phyglm <- function(formula, data, phy,
+                          Vx=NULL, Vy = NULL, times=2,
+                          distrib="uniform", btol=50, track=TRUE,...){
   #Error check
   if(class(formula)!="formula") stop("formula must be class 'formula'")
   if(class(data)!="data.frame") stop("data must be class 'data.frame'")
   if(class(phy)!="phylo") stop("phy must be class 'phylo'")
   if(distrib=="normal") warning ("distrib=normal: make sure that standard deviation 
                                  is provided for Vx or Vy")
-  else
-
   #Matching tree and phylogeny using utils.R
   datphy<-match_dataphy(formula,data,phy)
   full.data<-datphy[[1]]
@@ -77,8 +72,6 @@ intra_phyloglm <- function(formula,data,phy,
   if(!is.null(Vx) && sum(is.na(full.data[,Vx]))!=0){
     full.data[is.na(full.data[,Vx]),] <- 0}
 
-
-
   #Function to pick a random value in the interval
   if (distrib=="normal") funr <- function(a, b) {stats::rnorm(1,a,b)}
   else  funr <- function(a,b) {stats::runif(1,a-b,a+b)}
@@ -133,8 +126,8 @@ intra_phyloglm <- function(formula,data,phy,
       optpar               <- mod$alpha
 
 
-      if(track==TRUE) print(paste("intra: ",i,sep=""))
-      
+      if(track==TRUE) cat("\r","Simu = ", i," / ", times)
+
       #write in a table
       estim.simu <- data.frame(i, intercept, se.intercept, pval.intercept,
                                slope, se.slope, pval.slope, aic.mod, optpar,

diff --git a/R/intra_phylolm.R b/R/intra_phylolm.R
@@ -45,93 +45,79 @@
 #' @return \code{stats}: Statistics for model parameters. \code{sd_intra} is the standard deviation 
 #' due to intraspecific variation. \code{CI_low} and \code{CI_high} are the lower and upper limits 
 #' of the 95% confidence interval.
-#' @examples
-#' \dontrun{
-#' library(sensiPhy)
-#' 
-#' # Loading data and phylogeny:
-#' data(alien)
-#' trait <- log10(alien$data[,-1]+1)
-#' phy <- alien$phy[[1]]
-#' 
-#' # Running 50 models generating random predictor values 
-#' with a normal distribution  
-#' mods<-intra_phylm(mass~gesta,trait,phy,Vx="SD_gesta",times=50)
-#' summary(mods)
-#' sensi_plot(mods)
-#' }
-#' 
 #' @author Caterina Penone & Pablo Ariel Martinez
 #' @seealso \code{\link{sensi_plot}}
 #' @references Here still: reference to phylolm paper + our own?
 #' @export
 
 
-intra_phylm <- function(formula,data,phy,
-                        Vy=NULL,Vx=NULL,
-                        times=2,distrib="normal",
-                        model="lambda",track=TRUE,...){
+intra_phylm <- function(formula, data, phy,
+                        Vy = NULL, Vx = NULL,
+                        times = 2, distrib = "normal",
+                        model = "lambda", track = TRUE, ...){
   #Error check
-  if(class(formula)!="formula") stop("formula must be class 'formula'")
-  if(class(data)!="data.frame") stop("data must be class 'data.frame'")
-  if(class(phy)!="phylo") stop("phy must be class 'phylo'")
-  if(distrib=="normal") warning ("distrib=normal: make sure that standard deviation 
+  if(class(formula) != "formula") stop("formula must be class 'formula'")
+  if(class(data) != "data.frame") stop("data must be class 'data.frame'")
+  if(class(phy) != "phylo") stop("phy must be class 'phylo'")
+  if(distrib == "normal") warning ("distrib=normal: make sure that standard deviation 
                                  is provided for Vx or Vy")
-  else
-
+
   #Matching tree and phylogeny using utils.R
-  datphy<-match_dataphy(formula,data,phy)
-  full.data<-datphy[[1]]
-  phy<-datphy[[2]]
+  datphy <- match_dataphy(formula, data, phy)
+  full.data <- datphy[[1]]
+  phy <- datphy[[2]]
 
-  resp<-all.vars(formula)[1]
-  pred<-all.vars(formula)[2]
+  resp <- all.vars(formula)[1]
+  pred <- all.vars(formula)[2]
 
-  if(!is.null(Vy) && sum(is.na(full.data[,Vy]))!=0){
-    full.data[is.na(full.data[,Vy]),] <- 0}
+  if(!is.null(Vy) && sum(is.na(full.data[, Vy])) != 0) {
+    full.data[is.na(full.data[, Vy]), ] <- 0}
 
-  if(!is.null(Vx) && sum(is.na(full.data[,Vx]))!=0){
-    full.data[is.na(full.data[,Vx]),] <- 0}
+  if(!is.null(Vx) && sum(is.na(full.data[, Vx])) != 0) {
+    full.data[is.na(full.data[, Vx]), ] <- 0}
 
 
 
   #Function to pick a random value in the interval
-  if (distrib=="normal") funr <- function(a,b) {stats::rnorm(1,a,b)}
-  else  funr <- function(a,b) {stats::runif(1,a-b,a+b)}
+  if (distrib == "normal") funr <- function(a,b) {stats::rnorm(1,a,b)}
+  else  funr <- function(a,b) {stats::runif(1, a - b, a + b)}
 
   #Create the results data.frame
-  intra.model.estimates<-data.frame("n.intra"=numeric(),"intercept"=numeric(),"se.intercept"=numeric(),
-                                    "pval.intercept"=numeric(),"slope"=numeric(),"se.slope"=numeric(),
-                                    "pval.slope"=numeric(),"aic"=numeric(),"optpar"=numeric())
-
-
+  intra.model.estimates <- data.frame("n.intra" = numeric(),"intercept" = numeric(),
+                                    "se.intercept" = numeric(), 
+                                    "pval.intercept" = numeric(),
+                                    "slope" = numeric(),
+                                    "se.slope" = numeric(),
+                                    "pval.slope" = numeric(),"aic" = numeric(),
+                                    "optpar" = numeric())
   #Model calculation
-  counter=1
+  counter = 1
   errors <- NULL
-  c.data<-list()
+  c.data <- list()
 
   for (i in 1:times) {
-
     ##Set response and predictor variables
     #Vy is not provided or is not numeric, do not pick random value
-    if(!inherits(full.data[,resp], c("numeric","integer")) || is.null(Vy)) {full.data$respV<-full.data[,resp]}
+    if(!inherits(full.data[,resp], c("numeric","integer")) || is.null(Vy)) 
+        {full.data$respV <- model.frame(formula, data = full.data)[,1]}
 
     #choose a random value in [mean-se,mean+se] if Vy is provided
-    if(!is.null(Vy))
-    {full.data$respV<-apply(full.data[,c(resp,Vy)],1,function(x)funr(x[1],x[2]))}
-
+    if (!is.null(Vy))
+    {full.data$respV <- apply(full.data[,c(resp,Vy)],1,function(x)funr(x[1],x[2]))}
 
     #Vx is not provided or is not numeric, do not pick random value
-    if(!inherits(full.data[,pred], c("numeric","integer")) || is.null(Vx)) {full.data$predV<-full.data[,pred]}
+    if (!inherits(full.data[,pred], c("numeric","integer")) || is.null(Vx))
+        {full.data$predV <- model.frame(formula, data = full.data)[,2]}
 
     #choose a random value in [mean-se,mean+se] if Vx is provided
     if(!is.null(Vx))
-    {full.data$predV<-apply(full.data[,c(pred,Vx)],1,function(x)funr(x[1],x[2]))}
+    {full.data$predV <- apply(full.data[,c(pred,Vx)],1,function(x)funr(x[1],x[2]))}
 
     #model
-    mod = try(phylolm::phylolm(respV~predV, data=full.data, model=model,phy=phy),FALSE)
+    mod = try(phylolm::phylolm(respV ~ predV, data = full.data, model = model,
+                               phy = phy), FALSE)
 
-    if(isTRUE(class(mod)=="try-error")) {
+    if(isTRUE(class(mod) == "try-error")) {
       error <- i
       names(error) <- rownames(c.data$full.data)[i]
       errors <- c(errors,error)
@@ -147,17 +133,14 @@ intra_phylm <- function(formula,data,phy,
       pval.slope           <- phylolm::summary.phylolm(mod)$coefficients[[2,4]]
       aic.mod              <- mod$aic
       n                    <- mod$n
-      #d                    <- mod$d
 
-      if (model == "BM"){
+      if (model == "BM") {
         optpar <- NA
       }
-      if (model != "BM"){
-        optpar               <- mod$optpar
+      if (model != "BM") {
+        optpar <- mod$optpar
       }
-
-      if(track==TRUE) print(paste("intra: ",i,sep=""))
-
+      if(track == TRUE) cat("\r","Simu = ", i," / ", times)
       #write in a table
       estim.simu <- data.frame(i, intercept, se.intercept, pval.intercept,
                                slope, se.slope, pval.slope, aic.mod, optpar,
@@ -170,20 +153,21 @@ intra_phylm <- function(formula,data,phy,
 
   #calculate mean and sd for each parameter
   #variation due to intraspecific variability
-  mean_by_randomval<-stats::aggregate(.~n.intra, data=intra.model.estimates, mean)
+  mean_by_randomval <- stats::aggregate(.~n.intra, data = intra.model.estimates,
+                                        mean)
 
-  statresults<-data.frame(min=apply(intra.model.estimates,2,min),
-                          max=apply(intra.model.estimates,2,max),
-                          mean=apply(intra.model.estimates,2,mean),
-                          sd_intra=apply(mean_by_randomval,2,stats::sd))[-1,]
+  statresults <- data.frame(min = apply(intra.model.estimates, 2, min),
+                          max = apply(intra.model.estimates, 2, max),
+                          mean = apply(intra.model.estimates, 2, mean),
+                          sd_intra = apply(mean_by_randomval, 2, stats::sd))[-1, ]
 
   statresults$CI_low  <- statresults$mean - qt(0.975, df = times-1) * statresults$sd_intra / sqrt(times)
   statresults$CI_high <- statresults$mean + qt(0.975, df = times-1) * statresults$sd_intra / sqrt(times)
 
-  res <- list(formula=formula,
-              datas=full.data,
-              model_results=intra.model.estimates,N.obs=n,
-              stats=statresults)
+  res <- list(formula = formula,
+              datas = full.data,
+              model_results = intra.model.estimates, N.obs = n,
+              stats = statresults)
   class(res) <- "sensiIntra"
   return(res)
 }