mocasin: Function which computes selective p-values and intervals for...

In davidruegamer/selfmade: SELective inference For Mixed and ADditive model Estimators

Function which computes selective p-values and intervals for gamm4 and merMod objects

Description

Function which computes selective p-values and intervals for gamm4 and merMod objects

Usage

mocasin(mod, checkFun, this_y = NULL, nrSamples = 1000,
  bayesian = TRUE, varInTestvec = c("est", "minMod", "varY",
  "supplied"), varForSampling = c("est", "minMod", "varY", "supplied"),
  VCOV_vT = NULL, VCOV_sampling = NULL, conditional = TRUE,
  name = NULL, nrlocs = 7, which = NULL, vT = NULL, G = NULL,
  efficient = TRUE, trace = TRUE)

Arguments

mod	an object of class merMod or result of gamm4 function
checkFun	a function of y, a vector of the same length as the original response vector which returns TRUE or FALSE depending on whether the selection event for a given y corresponds to the original model selection. See the example for more details.
this_y	original response vector (explicit reference may be necessary for certain model classes)
nrSamples	integer; the number of Monte Carlo samples to be used for inference (defaults to 1000)
bayesian	logical; whether or not to use a bayesian type covariance
varInTestvec	for expert use only; variance used in the test vector definition
varForSampling	variance used for inference; per default the estimated variance of mod is used. Other options are a conservative estimate based on the variance of the response is used ("varY") or to supply a numeric value to base inference on a customize variance
VCOV_vT	for expert use only; VCOV used in the test vector definition
VCOV_sampling	covariance matrix of dimension of the response used for inference; per default the estimated covariance of mod is used. Otherwise a matrix must be supplied on which basis inference is conducted. If the true covariance is unknown, an conservative alternative to plugging in the estimator is given by using the covariance of the refitted mixed model, for which all fixed effects but the intercept are excluded.
conditional	logical; determines whether to use the conditional or marginal approach when mod is of class merMod, i.e., inference is sought for a linear mixed model
name	character; for the gamm4-case: the name of the covariate, for which inference is done
nrlocs	integer; for the gamm4-case: the number of locations tested for non-linear effects
which	integer; for the merMod-case: defining the effect for which inference is done
vT	list of vectors (optional); if inference is sought for a customized test vector, this argument can be used
G	true random effect covariance (optional)
efficient	logical; whether or not to compute the test statistic based on an (efficient) weighted LS estimator instead of a OLS estimator for the marginal model
trace	logical; if TRUE, a progress bar is printed in the console

Details

Note that the additive and conditional mixed model approach currently only works for a diagonal error covariance.

Examples

library(lme4)
if(require(lmerTest)){

##### BASED ON lmerTest HELP PAGE #########
# define function to fit a model based on response
modFun <- function(y)
{
  ham$y <- y
  lmer(y ~ Gender + Information * Product + (1 | Consumer) +
  (1 | Product), data=ham)
  
  }

# define a function to select a model (based on a model)
selFun <- function(mod) step(mod)

# define a function which extracts the results of the selection procedure
extractSelFun <- function(this_mod){

this_mod <- attr(this_mod, "model")
if(class(this_mod)=="lm") 
  return(attr(this_mod$coefficients, "names")) else
    return(c(names(fixef(this_mod)), names(getME(this_mod, "theta"))))

}


## backward elimination of non-significant effects:
(step_result <- selFun(modFun(ham$Informed.liking)))
attr(step_result, "model")
## Elimination tables for random- and fixed-effect terms:
(sel <- extractSelFun(step_result))

## Now we can finally define the function checking the congruency
## with the original selection

checkFun <- function(yb){ 

 this_mod <- modFun(yb)
 setequal( extractSelFun(selFun(this_mod)), sel )
 
 }

# Now let's compute valid p-values conditional on the selection
res <- mocasin(attr(step_result, "model"), this_y = ham$Informed.liking,
          checkFun = checkFun, which = 1:4, nrSamples = 50, trace = FALSE)

# print(res)
} 

# gamm4 example similar to the one from gamm4 help page
if(require(gamm4)){
set.seed(0) 
dat <- gamSim(1,n=500,scale=2) ## simulate 4 term additive truth

dat$y <- 3 + dat$x0^2 + rnorm(n=500)
br <- gamm4(y~ s(x0) + s(x1), data = dat)
summary(br$gam) ## summary of gam

# do not use any selection
checkFun <- function(yb) TRUE

res <- mocasin(br, this_y = dat$y,
          checkFun = checkFun, 
          nrlocs = c(0.7,1), 
          nrSamples = 100)
}

davidruegamer/selfmade documentation built on Feb. 2, 2023, 7:51 a.m.