debug

Author: livioivil

DESCRIPTION

@@ -8,4 +8,4 @@ Description: runs resampling-based tests jointly (e.i. sign-flip score tests [He
 Imports: flipscores, ggplot2, methods, flip 
 Encoding: UTF-8
 License: GPL-2
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2

NAMESPACE

@@ -5,6 +5,7 @@ S3method(p.adjust,jointest)
 S3method(plot,jointest)
 S3method(summary,jointest)
 export(combine)
+export(combine_contrasts)
 export(flip2sss)
 export(join_flipscores)
 import(dplyr)

R/combine.R

@@ -32,37 +32,37 @@ combine <- function (mods, comb_funct = "maxT", by = NULL, by_list=NULL, tail =
   res = lapply(1:length(combined), .npc2jointest, 
                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))
+  res=list(Tspace=jointest:::.get_all_Tspace(res),summary_table=jointest:::.get_all_summary_table(res))
   class(res) <- unique(c("jointest", class(res)))
   res
 }
 
 
 ###################
-#' (nonparametric) combination of jointest object 
-#' @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 comb_funct  Combining function to be used. Several functions are implemented: "mean", "median", "Fisher", "Liptak", (equal to) "Stoufer", "Tippet", (equal to) "minp", "maxT" (the default). 
-#' 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 tail direction of the alternative hypothesis. It can be "two.sided" (or 0, the default), "less" (or -1) or "greater" (or +1)
 #' @export
-# @describeIn combine \code{combine_contrasts} combines the tests derived from the contrasts of a factor variable to get a global test for the factor (i.e. categorical predictor). It has strong analogies with anova test.
+#' @describeIn combine \code{combine_contrasts} combines the tests derived from the contrasts of a factor variable to get a global test for the factor (i.e. categorical predictor). It has strong analogies with anova test.
 
 combine_contrasts <- function (mods, comb_funct = "maxT", tail = 0) 
 {
   # names(mods) = .set_mods_names(mods)
-  combined=paste(mods$summary_table$Model,mods$summary_table$.assign,sep="_")
-  if (is.null(combined)){
-    warning('There is not column $summary_table$.assign, nothing is done.' )
-  }
-  uniq_nm = unique(combined)
-  combined = lapply(uniq_nm, function(nm) which(combined == nm))
-  names(combined) = uniq_nm
-
-  res = lapply(1:length(combined), .npc2jointest, 
-               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)))
+  # combined=paste(mods$summary_table$Model,mods$summary_table$.assign,sep="_")
+  # if (is.null(combined)){
+  #   warning('There is not column $summary_table$.assign, nothing is done.' )
+  # }
+  # uniq_nm = unique(combined)
+  # combined = lapply(uniq_nm, function(nm) which(combined == nm))
+  # names(combined) = uniq_nm
+  # 
+  # res = lapply(1:length(combined), .npc2jointest, 
+  #              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
+  smr=apply(res$summary_table[,c("Model",".assign"),drop=FALSE],1,paste,collapse=".")
+  new_names=sapply(unique(smr),function(x).find_common_pattern(res$summary_table$Coeff[smr==x]))
+  res=combine(res,by=c("Model",".assign"),comb_funct = comb_funct, tail = tail)
+  res$summary_table$Coeff=new_names
   res
 }
 
@@ -72,15 +72,15 @@ combine_contrasts <- function (mods, comb_funct = "maxT", tail = 0)
   comb_name = names(combined)[id]
   if (is.null(comb_name)) 
     comb_name = "Combined"
-  Tspace = npc(mods$Tspace[, combined[[id]],drop=FALSE], comb_funct = comb_funct, 
+  Tspace = jointest:::npc(mods$Tspace[, combined[[id]],drop=FALSE], comb_funct = comb_funct, 
                tail = tail)
   colnames(Tspace) = comb_name
   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))
+                             p = jointest:::.t2p_only_first(Tspace, tail = 1))
   list(Tspace = Tspace, summary_table = summary_table)
 }
 

R/flip2sss.R

@@ -5,30 +5,32 @@
 #' @param formula A formula or a list of formulas. It can be a complete model as.formula or a list of formulas, one for each element produced by the function.
 #' @param data The dataset to be used for fitting the model.
 #' @param cluster A vector or a formula evaluated on the data that defines the clusters.
+#' @param family as in \code{glm}, but given as a character. Not used if argument \code{summstats_within} is not \code{NULL}.
 #' @param summstats_within A vector of summary statistics model within the data or a function with argument data.
 #' @param ... Other arguments passed to the `flipscores` function.
-#'
+#' @export
 #' @return A jointest object containing the model results. Note that the flipscores models for each coefficient within are also included in the jointest object.
 #'
 #' @examples
-#' set.seed(1)
-#' N=50
-#' n=rpois(N,10)
+#' set.seed(123)
+#' N=20
+#' n=rpois(N,20)
 #' reff=rep(rnorm(N),n)
 #' 
 #' D=data.frame(X1=rnorm(length(reff)),
 #'              X2=rep(rnorm(N),n),
 #'              Grp=factor(rep(rep(LETTERS[1:3],length.out=N),n)),
 #'              SOGG=rep(1:N,n))
-#' D$Y=rbinom(n=nrow(D),prob=dlogis((D$Grp=="B") * D$X1 + D$X2),size=1)
+#' D$Y=rbinom(n=nrow(D),prob=plogis( 2*D$X1 * (D$Grp=="B") +  2*D$X2+reff),size=1)
 #' 
 #' formula <- Y ~ Grp * X1 + X2
 #' cluster <- factor(D$SOGG)
 #' library(logistf)
 #' summstats_within <- 'logistf::logistf(Y ~ X1, family = binomial(link = "logit"),control=logistf::logistf.control(maxit=100))'
 #' #summstats_within <- 'glm(Y ~ X1, family = binomial(link = "logit"))'
 #' library(jointest)
-#' res <- flip2sss(formula, D, cluster, summstats_within)
+#' res <- flip2sss(formula, D, cluster, summstats_within=summstats_within)
+#' res <- flip2sss(formula, D, cluster, family="binomial")
 #' summary(res)
 #' summary(combine(res))
 #' summary(combine(res,by="Model"))
@@ -40,6 +42,7 @@
 flip2sss <- function(formula=NULL,
                      data=NULL,
                      cluster=NULL,
+                     family="gaussian",
                      summstats_within=NULL,
                      ...){
 
@@ -50,11 +53,11 @@ flip2sss <- function(formula=NULL,
 
   ###################
 
-  between_objs = .get_sets_vars_between(formula, data, cluster)  
+  vars_between_within = .get_vars_between_within(formula, data, cluster)  
 
-  vars_between=between_objs$vars_between
-  data2lev=between_objs$data2lev
-  rm(between_objs)
+  vars_between=vars_between_within$vars_between
+  vars_within=vars_between_within$vars_within
+  rm(vars_between_within)
   ## make the second level dataset  
   set_between = unique(unlist(vars_between))
   vars_between_formulas = lapply(vars_between, paste0, collapse = "+")
@@ -63,13 +66,15 @@ flip2sss <- function(formula=NULL,
   names(vars_between_formulas) = names(vars_between)
 
   #####################
+  if(is.null(summstats_within))
+    summstats_within=paste0("glm(",formula[[2]],formula[[1]],vars_within,", family=",family,")")
+  set_between=c(".cluster",set_between[set_between!="1"])
   data$.cluster=cluster
-  ss_within = data %>% 
-    group_by(.cluster) %>%
+  data2lev = data %>% 
+    group_by(data[set_between]) %>%
     summarise(as.data.frame(t(coefficients(eval(parse(text=summstats_within))))))
-  ss_within$.cluster=NULL
-  names(ss_within) = gsub("\\W", ".", names(ss_within))
-  data2lev=cbind(data2lev,ss_within)
+  data2lev$.cluster=NULL
+  names(data2lev) = gsub("\\W", ".", names(data2lev))
 
   mods = lapply(vars_between_formulas, function(frm) glm(eval(frm, parent.frame()), data = data2lev))
 
@@ -92,44 +97,62 @@ flip2sss <- function(formula=NULL,
 }
 
 ######################
-.get_sets_vars_between <- function(formula, data, cluster){
+.expand_form <- function(FUN){
+  out <- reformulate(labels(terms(FUN)), FUN[[2]])
+  out
+}
+######################
+.get_vars_between_within <- function(formula, data, cluster){
+  formula=.expand_form(formula)
   #clst_vals = unique(cluster) NON SERVE!
   D = model.matrix(formula, data = data)
 
   ## find constant cols within cluster
   const_id = do.call(rbind, by(D, cluster, function(D) as.data.frame(t(apply(D, 2, is.constant)))))
   vars_between_intercept = apply(const_id, 2, all)
 
-  data2lev=by(D,cluster, function(d) d[1,vars_between_intercept,drop=FALSE])
-  data2lev=do.call(rbind,data2lev)
-  
-  terms = colnames(D)#attributes(terms(formula))$term.labels#labels(terms(reformulate(as.character(formula))))
-  ids_const = unique(attributes(D)$assign[vars_between_intercept])
-  # attributes(D)$assign%in%ids_const
-   if(attributes(terms(formula))$intercept==1){
-    ids_const = ids_const + 1
-    terms = c("1", terms)
-   }
-  intercept_vars_between_names=names(vars_between_intercept[vars_between_intercept])
-  intercept_vars_between_names=gsub("^\\(Intercept\\)$","1",intercept_vars_between_names)
-  vars_between_names = list(.Intercept. = intercept_vars_between_names)
-  
-  cors = suppressWarnings(by(D[, !vars_between_intercept], cluster, cor))
-  dim3 = c(dim(cors[[1]]), length(cors))
-  cors = array(unlist(cors), dim3)
-  cors = suppressWarnings(apply(cors, c(1, 2), min, na.rm = TRUE))
-  cors[is.infinite(cors)]=1
-  ei = eigen(cors)
-  vars_between_others = lapply(which(ei$values > .1), function(cls) names(vars_between_intercept[!vars_between_intercept])[which(ei$vectors[, cls] != 0)])
-  within_coeffs = sapply(vars_between_others, function(preds_vars) Reduce(intersect, strsplit(preds_vars, ":")))  
-  
-  vars_between_others = lapply(1:length(within_coeffs), function(i) gsub(within_coeffs, "1", vars_between_others[[i]]))
-  vars_between_others = lapply(1:length(within_coeffs), function(i) gsub(":1$", "", vars_between_others[[i]]))
-  vars_between_others = lapply(1:length(within_coeffs), function(i) gsub("^1:", "", vars_between_others[[i]]))
-  
-  names(vars_between_others) = within_coeffs
-  vars_between_names = c(vars_between_names, vars_between_others)
-  list(data2lev=data2lev,vars_between_names=vars_between_names)
-}
+  #vars_between_intercept
+  between_vars=attributes(D)$assign[vars_between_intercept]
+  intercept=0%in%between_vars
+  between_vars=unique(between_vars)
+  between_vars_intercept<-
+    between_vars<-attr(terms(formula),"term.labels")[between_vars]
+  if(intercept)
+    between_vars_intercept=c("1",between_vars_intercept)
 
+  # formula_preds=paste(collapse = "+",between_vars)
+  # formula_preds_between_intercept=paste0(formula[[2]],formula[[1]],formula_preds)
+  
+  # data_intercept=data[,between_vars]
+  
+  
+  
+  ############# within variables
+  within_vars=attributes(D)$assign[!vars_between_intercept]
+  within_vars=unique(within_vars)
+  within_vars_all=attr(terms(formula),"term.labels")[within_vars]
+  within_vars=within_vars_all
+  for(x in between_vars)
+    within_vars=gsub(x,"",within_vars)
+  within_vars=gsub(":","",within_vars)
+  within_vars=unique(within_vars)
+  
+  ############# others between variables
+  vars_between=lapply(within_vars, function(x) {
+    temp=within_vars_all[grep(x,within_vars_all)]
+    temp=gsub(x,"",within_vars_all)
+    temp=gsub(":","",temp)
+    if(any(temp=="")) temp[temp==""]="1"
+    temp
+  })
+  names(vars_between)=within_vars
+  
+  temp=list(".Intercept."=between_vars_intercept)
+  vars_between=c(temp,vars_between)
+  
+  
+  list(vars_within=within_vars,
+       vars_between=vars_between)
+}
+ 
 is.constant <- function(x) length(unique(x)) == 1