mhtp: Multiple testing procedure for variable selection in linear...

In MMS: Fixed Effects Selection in Linear Mixed Models

Description

Performs a multiple hypotheses testing in linear mixed models

Usage

mhtp(data,Y,z,grp,D,fix,rand,alpha,step,num,ordre,m,show,IT,maxq,speed)

Arguments

data	Input matrix of dimension n * p; each row is an observation vector. The intercept should be included in the first column as (1,...,1). If not, it is added.
Y	Response variable of length n.
z	Random effects matrix. Of size n*q.
grp	Grouping variable of length n.
D	Logical value. If TRUE, the random effects are considered to be independent, i.e. Psi is a diagonal matrix. D=TRUE should be used with nested grouping factors.
fix	Number of variables which are not submitted to selection. They have to be in the first columns of data. Default is 1, the selection is not performed on the intercept.
rand	A vector of length q: each entry k is the position of the random effects number k in the data matrix, 0 otherwise. If z contains variables that have both a fixed and a random effect, it is advised to not submit them to selection.
alpha	A user supplied type I error sequence. Default is (0.1,0.05).
step	The algorithm performs at most step iterations. Default is 3000.
num	Number of variables one wishes to order. Default is min(n-1,p-1,30).
ordre	Several possible algorithms to order the variables, ordre=c("bolasso","pval","pval_hd","FR"). "bolasso" uses the dyadic algorithm with the Bolasso techniaue, "pval" uses the p-values obtained with a regression on the full set of variables (only when p<n), "pval_hd" uses marginal regression, "FR" uses Forward Regression. Default is "bolasso".
m	Number of bootstrapped iteration of the Lasso. Only use if the algorithm is set to "bolasso". Default is m=10.
show	Vector of logical values, show=(showordre,showresult,showit). Default is (0,0,0). If showordre=TRUE, show the ordered variables at each step of the algorithm. if showresult=TRUE, show the value of the statistics and the estimated quantile at each step of the procedure; if ordre=bolasso. if showit=TRUE, show the iterations of the algorithm.
IT	Number of simulations in the calculation of the quantile. Default is 10000.
maxq	Number of maximum multiple hypotheses testing to do. Default is min(log(min(n,p)-1,2),5).
speed	Logical value. If TRUE, the algorithm is speeded up once the criteria convergence in beta and u are fulfilled.

Details

mhtp performs fixed effects selection in linear mixed models. It is a combination of the mht function from the mht-package with the algorithm used in the lassop function; so you might want to have a look at the help of mht.

Value

A 'mhtp' object is returned for which refit is available.

data	List of the user-data: the scaled matrix used in the algorithm, the first column being (1,...,1); Y; z and grp.
beta	The estimated vector of fixed effects coefficients. Each row concern a specific user level alpha.
fitted.values	Fitted values calculated with the fixed effects and the random effects.
u	Matrix with #alpha columns. Each column is the concatenation of the estimated random effects (u_1',...,u_q')' for the user level alpha.
Psi	Variance of the random effects. Matrix of dimension q*q.
sigma_e	Variance of the residuals.
it	Number of iterations of the algorithm.
quantile	Array of all the estimated quantiles calculated during the procedure.
ordrebeta	All different order that has been used during the procedure.
converge	Did the algorithm converge?
call	The call that produced this object.
arg	List of all the arguments of the function (used to refit the function).

See Also

mht, refit.mhtp

Examples

## Not run: 
N <- 20           # number of groups
p <- 20            # number of covariates (including intercept)
q <- 2            # number of random effect covariates
ni <- rep(6,N)    # observations per group
n <- sum(ni)   # total number of observations

grp <- factor(rep(1:N,ni)) # grouping variable
grp=rbind(grp,grp)

beta <- c(1,2,4,3,rep(0,p-3)) # fixed-effects coefficients
x <- cbind(1,matrix(rnorm(n*p),nrow=n)) # design matrix

u1=rnorm(N,0,sd=sqrt(2))
u2=rnorm(N,0,sd=sqrt(2))
bi1 <- rep(u1,ni) 
bi2 <- rep(u2,ni)
bi <- rbind(bi1,bi2)

z=x[,1:2,drop=FALSE]
   
epsilon=rnorm(120)
y <- numeric(n)
for (k in 1:n) y[k] <- x[k,]%*%beta + t(z[k,])%*%bi[,k] + epsilon[k]

########
fit=mhtp(x,y,z,grp,D=0,fix=1,rand=c(1,2),alpha=0.1,num=15)
#fit=mhtp(x,y,z,grp,D=0,fix=1,rand=c(1,2),alpha=0.1,num=15,show=c(1,1,1))

## End(Not run)

MMS documentation built on May 2, 2019, 12:38 p.m.