...

Author: livio

.Rbuildignore

@@ -1 +1,3 @@
-(^\.github/.+$)|(^README.+$)
+(^\.github/.+$)|(^README.+$)|(^\docs/.+$)
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

NAMESPACE

@@ -1,10 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(as,jointest)
+S3method(p.adjust,jointest)
 S3method(plot,jointest)
 S3method(summary,jointest)
 export(combine)
 export(join_flipscores)
-export(p.adjust.fwer)
 import(flipscores)
 import(ggplot2)
 importFrom(flip,flip.adjust)

R/combine.R

@@ -4,52 +4,54 @@
 #' Alternativelly it can be a custom function that has a matrix as input. The function return a vector of length equal to the number of rows of the input matrix. 
 #' @param combined a) if \code{NULL} it combines all coefficients with the same names along models, 
 #' b) if a vector of characters, it gives a combination for each stats with these names 
-#' c) if a list of arrays, each array has two columns reporting model (first column) 
-#' and name of the coefficient to be combined.
+#' @param by vector of characters referring to the column's name in summary_table. Equal to "Coeff" by default.
 #' @param tail direction of the alternative hypothesis. It can be "two.sided" (or 0, the default), "less" (or -1) or "greater" (or +1)
 #' @export
 #' 
-combine <- function (mods, comb_funct = "maxT", combined = NULL, tail = 0) 
+combine <- function (mods, comb_funct = "maxT", combined = NULL, by=NULL, tail = 0) 
 {
-  names(mods) = .set_mods_names(mods)
-  smr = .get_all_tested_coeffs_names(mods)
+  # names(mods) = .set_mods_names(mods)
+  if(is.null(by)) by_cols="Coeff" 
+  smr=apply(res$summary_table[,by_cols,drop=FALSE],1,paste,collapse="_")    
   if (is.null(combined)) {
     uniq_nm = unique(smr)
     if (length(uniq_nm) == length(smr)) {
-      combined = list(Overall = smr)
+      combined = list(Overall = 1:ncol(res$Tspace))
     }
     else {
-      combined = lapply(uniq_nm, function(nm) which(smr == 
-                                                      nm))
+      combined = lapply(uniq_nm, function(nm) which(smr == nm))
       names(combined) = uniq_nm
     }
-  }
-  else if (!is.list(combined)) {
+  } else 
+    if (!is.list(combined)) {
     uniq_nm = combined
-    combined = lapply(uniq_nm, function(nm) which(smr == 
-                                                    nm))
+    combined = lapply(uniq_nm, function(nm) which(smr == nm))
     names(combined) = uniq_nm
   }
   else {
     uniq_nm = unique(unlist(combined))
   }
-  Tspace = .get_all_Tspace(mods)
+  
   res = lapply(1:length(combined), .npc2jointest, 
-               Tspace = Tspace, combined = combined, tail = tail, comb_funct = comb_funct)
+               mods = mods, combined = combined, tail = tail, comb_funct = comb_funct)
   names(res) = names(combined)
+  res=list(Tspace=.get_all_Tspace(res),summary_table=.get_all_summary_table(res))
   class(res) <- unique(c("jointest", class(res)))
   res
 }
 
-.npc2jointest <- function (id, Tspace, combined, tail, comb_funct) 
+.npc2jointest <- function (id, mods, combined, tail, comb_funct) 
 {
   comb_name = names(combined)[id]
   if (is.null(comb_name)) 
     comb_name = "Combined"
-  Tspace = npc(Tspace[, combined[[id]]], comb_funct = comb_funct, 
+  Tspace = npc(mods$Tspace[, combined[[id]]], comb_funct = comb_funct, 
                           tail = tail)
   colnames(Tspace) = comb_name
-  summary_table = data.frame(Coeff = comb_name, Stat = comb_funct, 
+  Coeff=mods$summary_table[combined[[id]],"Coeff"]
+  Coeff=unique(Coeff)
+  if(length(Coeff)>1) Coeff="many"
+  summary_table = data.frame(Coeff = Coeff, Stat = comb_funct, 
                              nMods = max(1, length(combined[[id]])), S = Tspace[1], 
                              p = .t2p_only_first(Tspace, tail = 1))
   list(Tspace = Tspace, summary_table = summary_table)

R/facilities.R

@@ -27,8 +27,9 @@
 #'       seed = 1, tested_coeffs = "X")
 #'summary(res)
 #'summary(combine(res))
-#'res=p.adjust.fwer(res)
+#'res=jointest:::p.adjust.jointest(res)
 #'summary(res)
+#'summary(combine(res,by="Model"))
 join_flipscores <- function (mods, tested_coeffs = NULL, n_flips = 5000, score_type = "standardized", 
                              statistics = "t",seed=NULL,...) 
 {
@@ -49,8 +50,10 @@ join_flipscores <- function (mods, tested_coeffs = NULL, n_flips = 5000, score_t
   n_obs=sapply(mods, function(mod) length(mod$y))
   n_obs=max(n_obs)
 
+  mods_names=names(mods)
+  
   FLIPS=flipscores:::.make_flips (n_obs=n_obs,n_flips=n_flips)
-    modflips = lapply(1:length(mods), function(i) {
+    mods = lapply(1:length(mods), function(i) {
       temp = flipscores(formula = eval(mods[[i]],parent.frame()), score_type = score_type, 
                         flips = eval(FLIPS), to_be_tested = tested_coeffs[[i]],
                         output_flips=FALSE
@@ -68,10 +71,14 @@ join_flipscores <- function (mods, tested_coeffs = NULL, n_flips = 5000, score_t
       temp$summary_table=.get_summary_table_from_flipscores(temp)
       temp
     })
-  if(is.null(names(mods))){
-    names(modflips)=paste0("mod",1:length(modflips))
+    
+  if(is.null(mods_names)){
+    names(mods)=paste0("mod",1:length(mos))
   } else
-    names(modflips) = names(mods)
-  class(modflips) <- unique(c("jointest", class(modflips)))
-  modflips
+    names(mods) = mods_names
+    
+
+  out=list(Tspace=.get_all_Tspace(mods),summary_table=.get_all_summary_table(mods))
+  class(out) <- unique(c("jointest", class(out)))
+  out
 }

R/join_flipscores.R

@@ -6,6 +6,7 @@
 # @importFrom MASS glm.nb
 # @importFrom plyr laply
 # @importFrom methods is
+# @importFrom methods as
 # @importFrom stats D as.formula model.matrix sd summary.glm update
 #' @examples
 #' n=20

R/jointest_package.R

@@ -57,24 +57,7 @@
 
 NULL
 # 
-# #' as.jointest method for a jointest object.
-# #' @rdname jointest-method
-# #' @param object an object of class \code{jointest}.
-# #' @param names_obj a vector of names, its length must be equal to the length of \code{object}
-# #' @param ... additional arguments to be passed
-# #' @method as jointest
-# #' @docType methods
-# #' @export
- 
-# as.jointest <- function (object, names_obj=NULL, ...) 
-# {
-#  TODO: calcolare summary_table in ogni elemento di object. se flipscores usa 
-#       .get_summary_table_from_flipscores()
-#   if(!is.null(names_obj)) names(object)=names_obj
-#   if (is.null(names(object))) names(object)=paste0("mod",1:length(object))
-#   class(object) <- unique(c("jointest", class(object)))
-#   object
-# }
+
 # #' print.jointest print method for a jointest object.
 # #' @param x a jointest object
 # #' @method print jointest
@@ -94,25 +77,37 @@ NULL
 
 summary.jointest <- function (object, ...) 
 {
-  if (length(object) > 1) {
-    summary_table = lapply(object, function(cmb) cmb$summary_table)
-    summary_table = do.call(rbind,summary_table)
-  }
-  else {
-    summary_table = object[[1]]$summary_table
-  }
-  rownames(summary_table) = names(object)
-  summary_table
+object$summary_table
 }
 
 .get_summary_table_from_flipscores <- function(object){
   tab = as.data.frame(summary(object)$coefficients)
   tab = tab[!is.na(tab[, "Score"]), ]
+  colnames(tab)[ncol(tab)]="p"
   tab = cbind( Coeff = rownames(tab), tab)
 }
 
 
 is_signif=NULL
+###########################
+#' as.jointest method for a jointest object.
+#' @rdname jointest-method
+#' @param object an object of class \code{jointest}.
+#' @param names_obj a vector of names, its length must be equal to the length of \code{object}
+#' @param ... additional arguments to be passed
+#' @method as jointest
+#' @docType methods
+#' @export
+
+as.jointest <- function (object, names_obj=NULL, ...)
+{
+  # TODO: calcolare summary_table in ogni elemento di object. se flipscores usa
+  #     .get_summary_table_from_flipscores()
+  if(!is.null(names_obj)) names(object)=names_obj
+  if (is.null(names(object))) names(object)=paste0("mod",1:length(object))
+  class(object) <- unique(c("jointest", class(object)))
+  object
+}
 #############################################
 #' plot.jointest summary method for a jointest object.
 #' @rdname jointest-method

R/methods.R

@@ -1,28 +1,29 @@
-#'Correct p-values of a jointest object for multiple testing (FamilyWise Error Rate)
-#'@param mods object of class \code{jointest} (it can be a list of glm or flipscores converted in a \code{jointest} object using \code{as.jointest})
-#'@param method any method implemented in \code{flip::flip.adjust} or a custom function. In the last case it must be a function that uses a matrix as input and returns a vector of adjusted p-values equal to the numbero fo columns of the inputed matrix.
-#'@param tail = 0
-#'@param ... further parameters
-#'@export
-#'@importFrom flip flip.adjust
-
-
-p.adjust.fwer <- function (mods, method = "maxT", tail = 0, ...) 
+# #'Correct p-values of a jointest object for multiple testing (i.e. Selective Inference)
+# #'@param mods object of class \code{jointest} (it can be a list of glm or flipscores converted in a \code{jointest} object using \code{as.jointest})
+# #'@param method any method implemented in \code{flip::flip.adjust} or a custom function. In the last case it must be a function that uses a matrix as input and returns a vector of adjusted p-values equal to the numbero fo columns of the inputed matrix.
+# #'@param tail = 0
+# #'@param ... further parameters
+# #'@export
+# #'@importFrom flip flip.adjust
+#' 
+#' p.adjust.jointest method for a jointest object.
+#' @rdname jointest-method
+#' @param object an object of class \code{jointest}.
+#' @param method any method implemented in \code{flip::flip.adjust} or a custom function. In the last case it must be a function that uses a matrix as input and returns a vector of adjusted p-values equal to the numbero fo columns of the inputed matrix.
+#' @param ... additional arguments to be passed
+#' @method  p.adjust jointest
+#' @docType methods
+#' @importFrom flip flip.adjust
+#' @export
+#' 
+p.adjust.jointest <- function (mods, method = "maxT", tail = 0, ...) 
 {
-  Tspace = as.matrix(.get_all_Tspace(mods))
-  colnames(Tspace) = paste0("v", 1:ncol(Tspace))
   if(is.character(method)){
-    p.adj = flip.adjust(.set_tail(Tspace, tail = tail), 
+    p.adj = flip.adjust(.set_tail(mods$Tspace, tail = tail), 
                             method = method) } 
   else if(is.function(method)){
-    p.adj = method(.set_tail(Tspace, tail = tail))
-                              
+    p.adj = method(.set_tail(mods$Tspace, tail = tail))
   }
-  tmp = lapply(1:length(mods), function(id) {
-    mods[[id]]$summary_table$p.adj = p.adj[id]
-    mods[[id]]
-  })
-  names(tmp)=names(mods)
-  class(tmp) <- unique(c("jointest", class(tmp)))
-  tmp
+  mods$summary_table$p.adj = p.adj
+  mods
 }