join_flipscores {jointest} | R Documentation |
join tests from multiverse models
Description
The function allows hypothesis testing across all plausible multiverse models ensuring strong family-wise error rate control.
Usage
join_flipscores(mods, tested_coeffs = NULL, n_flips = 5000,
score_type = "standardized", statistics = "t", seed=NULL, output_models =TRUE, ...)
Arguments
mods |
list of |
tested_coeffs |
list of the same length of |
n_flips |
number of flips, default 5000 |
score_type |
any valid type for |
statistics |
"t" is the only method implemented (yet). Any other value will not modify the score (a different statistic will only affect the multivariate inference, not the univariate one). |
seed |
|
output_models |
|
... |
any other further parameter. |
Value
A jointest
object, i.e., a list containing the following objects:
- Tspace
data.frame
where rows represents the sign-flipping transformed (plus the identity one) test and columns the variables.- summary_table
data.frame
containing for each model the estimated parameter(s), score(s), std error(s), test(s), partial correlation(s) and p-value(s).- mods
List of
glm
s orflipscores
objects.
Examples
library(jointest)
set.seed(123)
#EXAMPLE 1: Simulate data:
n=20
D=data.frame(X=rnorm(n),Z1=rnorm(n),Z2=rnorm(n))
D$Y=D$Z1+D$X+rnorm(n)
# Run four glms abd combine it in a list
mod1=glm(Y~X+Z1+Z2,data=D)
mod2=glm(Y~X+poly(Z1,2)+Z2,data=D)
mod3=glm(Y~X+poly(Z1,2)+poly(Z2,2),data=D)
mod4=glm(Y~X+Z1+poly(Z2,2),data=D)
mods=list(mod1=mod1,mod2=mod2,mod3=mod3,mod4=mod4)
# flipscores jointly on all models
res=join_flipscores(mods,n_flips = 1000)
summary(combine(res))
summary(combine(res, by="Model"))
summary(combine_contrasts(res))
#Simulate multivariate (50) bionomial responses
set.seed(123)
n=30
D=data.frame(X=rnorm(n),Z=rnorm(n))
Y=replicate(50,rbinom(n,1,plogis(.5*D$Z+.5*D$X)))
colnames(Y)=paste0("Y",1:50)
D=cbind(D,Y)
mods=lapply(1:50,function(i)eval(parse(text=
paste(c("glm(formula(Y",i,"~X+Z),data=D,family='binomial')"),collapse=""))))
# flipscores jointly on all models
res=join_flipscores(mods,n_flips = 1000,tested_coeffs ="X")
summary(res)
res=p.adjust(res)
summary(res)
# Compute lower bound for the true discovery proportion. See packages pARI and sumSome
# install.packages("sumSome")
# install.packages("pARI")
# library(sumSome)
# library(pARI)
# pARI returns a lower bound equals 0.24, i.e., at least 24% of the models
# have a significant effect related to X
# pARI::pARI(ix = c(1:50),pvalues = t(jointest:::.t2p(res$Tspace)),family = "simes",delta = 9)$TDP
# sumSome returns a lower bound equals 0.42, i.e., at least 42% of the models
# have a significant effect related to X
# sumSome::tdp(sumSome::sumStats(G = as.matrix(res$Tspace)))