flip2sss function for second stage summary statistics via flipscores

Author: livioivil

NAMESPACE

@@ -5,9 +5,12 @@ S3method(p.adjust,jointest)
 S3method(plot,jointest)
 S3method(summary,jointest)
 export(combine)
+export(flip2sss)
 export(join_flipscores)
+import(dplyr)
 import(flipscores)
 import(ggplot2)
+import(magrittr)
 importFrom(flip,flip.adjust)
 importFrom(stats,median)
 importFrom(stats,qnorm)

R/anova_jointest.R

@@ -0,0 +1,91 @@
+#' anova.jointest
+#' @description This is the \code{anova} method for \code{jointest} object. Remark: it performs type III deviance decomposition as in \code{car::Anova}.
+#' @param object (the object) \code{glm} (or \code{jointest}) object with the model under the null hypothesis (i.e. the covariates, the nuisance parameters).
+#' @param model1 a \code{glm} (or \code{jointest}) or a \code{matrix} (or \code{vector}). If it is a \code{glm} object, it has the model under the alternative hypothesis. The variables in \code{model1} are the same variables in \code{object} plus one or more variables to be tested.  Alternatively, if
+#' \code{model1} is a \code{matrix}, it contains the tested variables column-wise.
+#' @param score_type The type of score that is computed. It can be "orthogonalized", "effective" or "basic".
+#' Default is "orthogonalized". "effective" and "orthogonalized" take into account the nuisance estimation. The default is \code{NULL}, in this case the value is taken from \code{object}.
+#' @param n_flips The number of random flips of the score contributions.
+#' When \code{n_flips} is equal or larger than the maximum number of possible flips (i.e. n^2), all possible flips are performed. 
+#' Default is 5000.
+#' @param id a \code{vector} identifying the clustered observations. If \code{NULL} (default) observations are assumed to be independent. NOTE: if \code{object} is a \code{jointest} and \code{model$flip_param_call$id} is not \code{NULL}, this is considered in the inference.
+#' @param ... other parameters allowed in \code{stats::anova}.
+#' @examples
+#' set.seed(1)
+#' 
+#'
+#'
+
+
+anova.jointest <- function(object,
+                             ...){
+  
+  if(is.null(object$x)||(length(object$x)==0)) object=update(object,x=TRUE)
+  
+  if(is.null(score_type)) score_type = object$score_type
+  
+  anova_temp=get("anova.glm", envir = asNamespace("stats"),
+                 inherits = FALSE)
+  
+  ## comparison of 2 nested models
+  if(!is.null(model1)){ 
+    scores=compute_scores(model0 = object,
+                          model1 = model1,
+                          score_type = score_type)
+    mf <- match.call(expand.dots = TRUE)
+    if(!is.null(mf$flip_param_call$id))
+      scores=rowsum(scores,group = id)
+    
+    Tspace=sapply(1:ncol(scores), function(id_col){
+      score1=scores[,id_col,drop=FALSE]
+      attributes(score1)$scale_objects=attributes(scores)$scale_objects[[id_col]]
+      attributes(score1)$score_type=attributes(scores)$score_type
+      as.matrix(.flip_test_no_pval(score1, precompute_flips = FALSE,
+                                   .score_fun = .score_std,n_flips = n_flips))
+    })
+    dst=mahalanobis_npc(Tspace)
+    
+    out_param=anova_temp(object,model1,test="Rao")
+    heading2=attributes(out_param)$heading[2]
+    out_param = out_param[-1,]
+    out_param = out_param[,-c(1:2,4)]
+    names(out_param)[2]="Score"
+    out_param[[2]]=dst[1]
+    names(out_param)[3]="Pr(>Score)"
+    out_param[[3]]=.t2p(dst)
+    rownames(out_param)[1]="Model 2 vs Model 1"
+    
+  } else   { ## one anova for all variables
+    
+    varlist <- attr(object$terms, "variables")
+    if (!is.matrix(object[["x"]])) 
+      object$x = model.matrix(object)
+    varseq <- attr(object$x, "assign")
+    
+    ############################
+    subsets_npc=lapply(unique(varseq[varseq!=0]),function(i)which(varseq==i))
+    
+    res=mahalanobis_npc_multi(ids_list = subsets_npc,permT = as.matrix(object$Tspace))
+    # flip::npc(ps@permT,comb.funct = "mahalanobist",subsets = subsets_npc)
+    # res@res[, 3]=res@res[, 3]*nrow(ps@permT)
+    
+    out_param = anova_temp(object,test="Rao")
+    heading2=paste0("Model: ",deparse1(formula(object)))
+    out_param = out_param[-1,]
+    out_param = out_param[,-(2:4)]
+    names(out_param)[2]="Score"
+    out_param[[2]]=res[1,]
+    names(out_param)[3]="Pr(>Score)"
+    out_param[[3]]=apply(res,2,.t2p)
+  } # closes if
+  #make up
+  title <- paste0("Analysis of Deviance Table (Type III test)", 
+                  "\n\nModel: ", 
+                  object$family$family, ", link: ", object$family$link,"\n") 
+  attr(out_param,"heading")=c(title,heading2)
+  
+  
+  attr(out_param,"heading")[[1]]=  paste(attr(out_param,"heading")[[1]],sep="",
+                                         "\nInference is provided by jointest approach (",object$flip_param_call$n_flips," sign flips).\n")
+  return(out_param)
+}

R/combine.R

@@ -45,7 +45,7 @@ combine <- function (mods, comb_funct = "maxT", combined = NULL, by=NULL, tail =
   comb_name = names(combined)[id]
   if (is.null(comb_name)) 
     comb_name = "Combined"
-  Tspace = npc(mods$Tspace[, combined[[id]]], comb_funct = comb_funct, 
+  Tspace = npc(mods$Tspace[, combined[[id]],drop=FALSE], comb_funct = comb_funct, 
                           tail = tail)
   colnames(Tspace) = comb_name
   Coeff=mods$summary_table[combined[[id]],"Coeff"]

R/desktop.ini

@@ -0,0 +1,108 @@
+#' flip2sss
+#'
+#' This function fits a model based on the provided formula and data, accounting for clusters and summary statistics within the model.
+#'
+#' @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 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.
+#'
+#' @return A jointest object containing the model results. Note that the flipscores models for each coefficient within are also included in the jointest object.
+#' @export
+#'
+#' @examples
+#' formula <- SALUTE ~ GENERE * TIME + PROV_BAMB + ETA_B_ARRIVO
+#' data <- adolong
+#' cluster <- adolong$SOGG
+#' library(logistf)
+#' temp=logistf::logistf(SALUTE ~ TIME, family = binomial(link = "logit"),data=data,control=logistf.control(maxit=100))
+#' summstats_within <- 'logistf::logistf(SALUTE ~ TIME, family = binomial(link = "logit"),control=logistf::logistf.control(maxit=100))'
+#' res <- flip2sss(formula, data, cluster, summstats_within)
+#'attr(model.matrix(mods[[1]]),"assign")
+#' @import dplyr
+#' @import magrittr
+#' 
+flip2sss <- function(formula=NULL,
+                     data=NULL,
+                     cluster=NULL,
+                     summstats_within=NULL,
+                     ...){
+  
+  # if(is(cluster,"formula")){
+  #   cluster=eval(cluster,data)
+  # }
+  
+ 
+  ###################
+  
+  vars_between = .get_sets_vars_between(formula, data, cluster)  
+  
+  ## make the second level dataset  
+  set_between = unique(unlist(vars_between))
+  vars_between_formulas = lapply(vars_between, paste0, collapse = "+")
+  vars_between_formulas = paste(names(vars_between_formulas), vars_between_formulas, sep = "~")
+  vars_between_formulas = as.list(vars_between_formulas)
+  names(vars_between_formulas) = names(vars_between)
+  
+  #####################
+  data2lev = data %>% 
+    group_by(adolong[,set_between[-1]]) %>%
+    summarise(as.data.frame(t(coefficients(eval(parse(text=summstats_within))))))
+  names(data2lev) = gsub("\\W", ".", names(data2lev))
+  
+  mods = lapply(vars_between_formulas, function(frm) glm(eval(frm, parent.frame()), data = data2lev))
+  
+  for(i in 1:length(mods)){
+    mods[[i]]$call$data = eval(data2lev)
+    mods[[i]]$call$formula = eval(as.formula(vars_between_formulas[[i]]))
+  }
+  
+  res = join_flipscores(mods, n_flips = 5000, seed = 1)
+  # summary(res)
+  
+  res$summary_table$Coeff = apply(res$summary_table[, 2:1], 1, paste, collapse = ":")
+  res$summary_table$Coeff = gsub(":\\.Intercept\\.$", "", res$summary_table$Coeff)
+  res$summary_table$Coeff = gsub("\\(Intercept\\):", "", res$summary_table$Coeff)
+  
+  # res$summary_table$Model = "flip2sss"
+  colnames(res$Tspace) = apply(res$summary_table[, 2:1], 1, paste, collapse = "_model.")
+  res$mods=mods
+  res
+}
+
+######################
+.get_sets_vars_between <- function(formula, data, cluster){
+  clst_vals = unique(cluster)
+  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)
+  
+  terms = labels(terms(reformulate(as.character(formula))))
+  ids_const = unique(attributes(D)$assign[vars_between_intercept])
+  if(0 %in% ids_const){
+    ids_const = ids_const + 1
+    terms = c("1", terms)
+  }
+  vars_between_names = list(.Intercept. = terms[ids_const])
+  
+  cors = suppressWarnings(by(D[, !vars_between_intercept], cluster, cor))
+  dim3 = c(dim(cors[[1]]), length(cors))
+  cors = array(unlist(cors), dim3)
+  cors = apply(cors, c(1, 2), min, na.rm = TRUE)
+  ei = eigen(cors)
+  vars_between_others = lapply(which(ei$values > .1), function(cls) terms[-ids_const][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)
+  vars_between_names
+}
+
+is.constant <- function(x) length(unique(x)) == 1

R/flip2sss.R

@@ -78,10 +78,36 @@ join_flipscores <- function (mods, tested_coeffs = NULL, n_flips = 5000, score_t
   } else
     names(mods) = mods_names
 
-
+    assign=sapply(mods,function(mod)attr(model.matrix(mod),"assign"))
+    vars_orig=sapply(mods,function(mod)terms(mod$formula[[3]]))
+    names(vars_orig)=mods_names
+    
+    for(i in 1:length(assign)){
+      if(min(assign[[i]])==0){
+        assign[[i]]=assign[[i]]+1
+        vars_orig[[i]]=c(".Intercept.",vars_orig[[i]])
+      }
+    }
+    
+    if(length(assign)>1){
+      for(i in 2:length(assign)){
+        cnt=max(assign[[i-1]])
+        assign[[i]]=cnt+assign[[i]]
+      }
+    } 
+    
+    
+    temp=lapply(1:length(vars_orig), 
+           function(i) paste0("mod_",names(vars_orig)[i],"_",vars_orig[[i]]))
+    names_vars_orig=unlist(temp)
+    assign=unlist(assign)
+    names(assign)=NULL
+    temp=lapply(unique(assign),function(i) which(assign==i))
+    names(temp)=names_vars_orig
   out=list(Tspace=.get_all_Tspace(mods),
            summary_table=.get_all_summary_table(mods),
            mods=mods)
+  attr(out$Tspace,"orig_var")=temp
   class(out) <- unique(c("jointest", class(out)))
   out
 }

R/join_flipscores.R

@@ -57,14 +57,18 @@
 
 NULL
 # 
+#' non funziona e non funziona neppure summary se lo uso: 
+#' #' print.jointest print method for a jointest object.
+#' #' @rdname jointest-method
+#' #' @param object an object of class \code{jointest}.
+#' #' @param ... additional arguments to be passed
+#' #' @method print jointest
+#' #' @docType methods
+#' #' @export
+#' print.jointest <- function (object, ...) {
+#'   #summary.jointest(object, ...)
+#' }
 
-# #' print.jointest print method for a jointest object.
-# #' @param x a jointest object
-# #' @method print jointest
-# #' @docType methods
-# #' @rdname jointest-method
-# #' @export
-# 
 # 
 
 #' summary.jointest summary method for a jointest object.
@@ -133,6 +137,7 @@ plot.jointest <- function(object,
   if(p.values =="raw"){
     if(!is.null(D$`Pr(>|z|)`)) D$p.vals=D$`Pr(>|z|)` else
       D$p.vals=D$`Pr(>|t|)`
+    if(is.null(D$p.vals)) D$p.vals=D$p
   } else
     if(p.values=="adjusted")
       D$p.vals=D$p.adj

R/methods.R

@@ -1,32 +1,40 @@
 
 npc <- function (Tspace, comb_funct = "Fisher", tail = 0) 
 {
-  if (!(comb_funct %in% c("mean", "median"))) 
-    Tspace = .set_tail(Tspace, tail = tail)
-  if (comb_funct %in% c("Fisher", "Liptak", "Stoufer", 
-                        "Tippet", "minp")) {
+  if(is.character(comb_funct)){
+    implemented_comb_funcs=c("mean", "median","Fisher", 
+                           "Liptak", "Stoufer", 
+                           "Tippet", "minp","Mahalanobis")
+    comb_funct=match.arg(comb_funct,implemented_comb_funcs)
+    if (comb_funct %in% c("mean", "median"))
+      Tspace = .set_tail(Tspace, tail = tail)
+    if (comb_funct %in% c("Fisher", 
+                            "Liptak", "Stoufer", 
+                            "Tippet", "minp"))
     Tspace = .t2p(Tspace, tail = 1)
-  }
-  if (comb_funct == "minp") 
-    Tspace = -Tspace
-  if (comb_funct == "mean") {
-    Tspace = rowMeans(Tspace)
-    Tspace = .set_tail(Tspace, tail = tail)
-  }
-  else if (comb_funct %in% c("minp", "maxT")) {
-    Tspace = .rowMax(Tspace)
-  }
-  else if (comb_funct == "median") {
-    Tspace = .rowMedians(Tspace)
-    Tspace = .set_tail(Tspace, tail = tail)
-  }
-  else if (comb_funct == "Fisher") {
-    Tspace = rowSums(.comb_funct_fisher(Tspace))
-  }
-  else if (comb_funct %in% c("Liptak", "Stoufer")) {
-    Tspace = rowSums(.comb_funct_liptak(Tspace))
-  }
-  else if (is.function(comb_funct)) {
+  
+    if (comb_funct == "minp") 
+      Tspace = -Tspace
+    if (comb_funct == "mean") {
+      Tspace = rowMeans(Tspace)
+      Tspace = .set_tail(Tspace, tail = tail)
+    }
+    else if (comb_funct %in% c("minp", "maxT")) {
+      Tspace = .rowMax(Tspace)
+    }
+    else if (comb_funct == "median") {
+      Tspace = .rowMedians(Tspace)
+      Tspace = .set_tail(Tspace, tail = tail)
+    }
+    else if (comb_funct == "Fisher") {
+      Tspace = rowSums(.comb_funct_fisher(Tspace))
+    }
+    else if (comb_funct %in% c("Liptak", "Stoufer")) {
+      Tspace = rowSums(.comb_funct_liptak(Tspace))
+    } else if (comb_funct %in% c("Mahalanobis")) {
+      Tspace = flipscores:::mahalanobis_npc(Tspace)
+    }
+  } else if (is.function(comb_funct)) {
     Tspace = comb_funct(Tspace)
     comb_funct = "custom"
   }

R/npc.R

@@ -5,12 +5,11 @@
   tab=lapply(mods,function(x){
     if(is.null(x$Tspace)){stop("At least one Tspace is missing")}
     x$Tspace})
+  nms = sapply(tab,colnames)
   tab=do.call(cbind,tab)
-  colnames(tab) = paste0("S", 1:ncol(tab))
-  # tab=mods[[1]]$Tspace
-  # for (i in 2:length(mods)){
-  #   tab=cbind(tab,mods[[i]]$Tspace)
-  # }
+  colnames(tab)=nms
+  if(is.null(colnames(tab))) colnames(tab)=paste0("S", 1:ncol(tab))
+  
    tab
 }