Group-Adaptive Elastic Net Penalised Generalised Linear Models

Description

Fit linear and logistic regression models penalised with group-adaptive elastic net penalties. The group penalties correspond to groups of covariates defined by a co-data group set. The method accommodates inclusion of unpenalised covariates and overlapping groups. See Van Nee et al. (2021) <arXiv:2101.03875>.

Details

See squeezy for example code.

Author(s)

Mirrelijn M. van Nee [aut, cre], Tim van de Brug [aut], Mark A. van de Wiel [aut]

Maintainer: Mirrelijn M. van Nee <m.vannee@amsterdamumc.nl>

References

Mirrelijn M. van Nee, Tim van de Brug, Mark A. van de Wiel, "Fast marginal likelihood estimation of penalties for group-adaptive elastic net.", 2021

Partial derivatives of -log(ML) of ridge penalised GLMs

Description

Returns the partial derivatives (w.r.t. 'loglambdas') of the minus log Laplace approximation (LA) of the marginal likelihood of ridge penalised generalised linear models. Note: currently only implemented for linear and logistic regression.

Usage

dminML.LA.ridgeGLM(loglambdas, XXblocks, Y, sigmasq = 1, 
                  Xunpen = NULL, intrcpt = TRUE, model, minlam = 0, 
                  opt.sigma = FALSE)

Arguments

Value

Partial derivatives of the Laplace approximation of the minus log marginal likelihood to the model parameters 'loglambdas';

For opt.sigma=FALSE: Gx1-dimensional vector for the G log(group ridge penalties).

For opt.sigma=TRUE (linear model only): (G+1)x1-dimensional vector for the partial derivative to log(sigmasq) (first element) and for the G log(group ridge penalties).

Examples

#Simulate toy data
n<-100
p<-300
X <- matrix(rnorm(n*p),n,p)
Y <- rnorm(n)
groupset <- list(1:(p/2),(p/2+1):p)
sigmahat <- 2
alpha <- 0.5
tauMR <- c(0.01,0.005)

XXblocks <- lapply(groupset, function(x)X[,x]%*%t(X[,x]))

#compute partial derivatives of the minus log marginal likelihood to the penalties only
dminML.LA.ridgeGLM(loglambdas = log(sigmahat/tauMR), 
                  XXblocks, Y, sigmasq = sigmahat, 
                  model="linear",opt.sigma=FALSE)

#additionally, compute the partial derivative to the linear regression noise parameter sigma^2
dminML.LA.ridgeGLM(loglambdas = log(c(sigmahat,sigmahat/tauMR)), 
                   XXblocks, Y, sigmasq = sigmahat, 
                   model="linear",opt.sigma=TRUE)

Marginal AIC of a multi-group, ridge penalised GLM

Description

Compute the marginal AIC for the marginal likelihood (ML) of multi-group, ridge penalised generalised linear models. Note: currently only implemented for linear and logistic regression.

Usage

mAIC.LA.ridgeGLM(loglambdas, XXblocks, Y, sigmasq = 1, 
                Xunpen = NULL, intrcpt = TRUE, model, minlam = 0)

Arguments

Value

Examples

#Simulate toy data
n<-100
p<-300
X <- matrix(rnorm(n*p),n,p)
Y <- rnorm(n)
groupset <- list(1:(p/2),(p/2+1):p)
sigmahat <- 2
alpha <- 0.5
tauMR <- c(0.01,0.005)

XXblocks <- lapply(groupset, function(x)X[,x]%*%t(X[,x]))

#compute the mAIC of a co-data model with multiple groups
mAIC.LA.ridgeGLM(loglambdas=log(sigmahat/tauMR), XXblocks=XXblocks, 
                 Y = Y, sigmasq = sigmahat, model="linear")

#compute the mAIC of a co-data agnostic model, i.e. only one group of covariates
mAIC.LA.ridgeGLM(loglambdas=log(sigmahat/median(tauMR)), 
                 XXblocks=list(X%*%t(X)), 
                 Y = Y, sigmasq = sigmahat, model="linear")

-log(ML) of ridge penalised GLMs

Description

Returns the Laplace approximation (LA) of the minus log marginal likelihood of ridge penalised generalised linear models. Note: currently only implemented for linear and logistic regression.

Usage

minML.LA.ridgeGLM(loglambdas, XXblocks, Y, sigmasq = 1, 
                  Xunpen = NULL, intrcpt = TRUE, model, minlam = 0, 
                  opt.sigma = FALSE)

Arguments

Value

Laplace approximation of the minus log marginal likelihood for the ridge penalised GLM with model parameters 'loglambdas' and 'sigmasq' (for linear regression).

Examples

#Simulate toy data
n<-100
p<-300
X <- matrix(rnorm(n*p),n,p)
Y <- rnorm(n)
groupset <- list(1:(p/2),(p/2+1):p)
sigmahat <- 2
alpha <- 0.5
tauMR <- c(0.01,0.005)

XXblocks <- lapply(groupset, function(x)X[,x]%*%t(X[,x]))

#compute minus log marginal likelihood
minML.LA.ridgeGLM(loglambdas = log(sigmahat/tauMR), 
                  XXblocks, Y, sigmasq = sigmahat, 
                  model="linear")

Visual posterior check of multivariate normality of the linear predictors

Description

Produce a qq-plot to visually check whether the assumption of multivariate normality of the linear predictors is valid for the data and model fit with 'squeezy'.

Usage

normalityCheckQQ(X,groupset,fit.squeezy,nSim=500)

Arguments

Value

The qqplot of the empirical versus theoretical quantiles is plotted. If ‘ggplot2’ is installed, the plot is returned as ‘ggplot’ object.

Examples

#Simulate toy data
n<-100
p<-300
X <- matrix(rnorm(n*p),n,p)
Y <- rnorm(n)
groupset <- list(1:(p/2),(p/2+1):p)
sigmahat <- 2
alpha <- 0.5
tauMR <- c(0.01,0.005)

#Fit group-regularised elastic net model with squeezy
fit.squeezy <- squeezy(Y,X,groupset,alpha=alpha,
                       lambdas=sigmahat/tauMR,sigmasq=sigmahat,
                       lambdaglobal=mean(sigmahat/tauMR))
                       
#Check qq-plot
normalityCheckQQ(X,groupset,fit.squeezy)

Fit a group-adaptive elastic net penalised linear or logistic model

Description

Estimate group-specific elastic net penalties and fit a linear or logistic regression model.

Usage

squeezy(Y, X, groupset, alpha = 1, model = NULL, X2 = NULL, 
        Y2 = NULL, unpen = NULL, intrcpt = TRUE, 
        method = c("ecpcEN", "MML", "MML.noDeriv", "CV"), 
        fold = 10, compareMR = TRUE, selectAIC = FALSE, fit.ecpc = NULL, 
        lambdas = NULL, lambdaglobal = NULL, lambdasinit = NULL, 
        sigmasq = NULL, ecpcinit = TRUE, SANN = FALSE, minlam = 10^-3, 
        standardise_Y = NULL, reCV = NULL, opt.sigma = NULL, 
        resultsAICboth = FALSE, silent=FALSE)

Arguments

Value

If ‘compareMR’=TRUE, multi-ridge model is returned as well:

If independent test set ‘X2’ is given, predictions and MSE are returned:

If ‘selectAIC’=TRUE, the multi-group or single-group model with best AIC is selected. Results in ‘betaApprox’, ‘a0Approx’, ‘lambdaApprox’ contain those results of the best model. Summary results of both models are included as well:

Author(s)

Mirrelijn M. van Nee, Tim van de Brug, Mark A. van de Wiel

References

Mirrelijn M. van Nee, Tim van de Brug, Mark A. van de Wiel, "Fast marginal likelihood estimation of penalties for group-adaptive elastic net", arXiv preprint, arXiv:2101.03875 (2021).

Examples

#####################
# Simulate toy data #
#####################
p<-100 #number of covariates
n<-50 #sample size training data set
n2<-100 #sample size test data set
G<- 5 #number of groups

taugrp <- rep(c(0.05,0.1,0.2,0.5,1),each=p/G) #ridge prior variance
groupIndex <- rep(1:G,each=p/G) #groups for co-data
groupset <- lapply(1:G,function(x){which(groupIndex==x)}) #group set with each element one group
sigmasq <- 2 #linear regression noise
lambda1 <- sqrt(taugrp/2) #corresponding lasso penalty
#A Laplace(0,b) variate can also be generated as the difference of two i.i.d.
#Exponential(1/b) random variables
betas <-   rexp(p, 1/lambda1) -  rexp(p, 1/lambda1) #regression coefficients
X <- matrix(rnorm(n*p),n,p) #simulate training data
Y <- rnorm(n,X%*%betas,sd=sqrt(sigmasq))
X2 <- matrix(rnorm(n*p),n,p) #simulate test data
Y2 <- rnorm(n,X2%*%betas,sd=sqrt(sigmasq))

###############
# Fit squeezy #
###############
#may be fit directly..
res.squeezy <- squeezy(Y,X,groupset=groupset,Y2=Y2,X2=X2,
                       model="linear",alpha=0.5)


  #..or with ecpc-fit as initialisation
  if(requireNamespace("ecpc")){
    res.ecpc <- ecpc::ecpc(Y,X, #observed data and response to train model
                     groupsets=list(groupset), #informative co-data group set
                     Y2=Y2,X2=X2, #test data
                     model="linear",
                     hypershrinkage="none",postselection = FALSE)
    res.squeezy <- squeezy(Y,X, #observed data and response to train model
                           groupset=groupset, #informative co-data group set
                           Y2=Y2,X2=X2, #test data
                           fit.ecpc = res.ecpc, #ecpc-fit for initial values
                           model="linear", #type of model for the response
                           alpha=0.5) #elastic net mixing parameter
  }



summary(res.squeezy$betaApprox) #estimated elastic net regression coefficients
summary(res.squeezy$betaMR) #estimated multi-ridge regression coefficients
res.squeezy$lambdaApprox #estimated group elastic net penalties
res.squeezy$tauMR #multi-ridge group variances
res.squeezy$MSEApprox #MSE group-elastic net model
res.squeezy$MSEMR #MSE group-ridge model

#once fit, quickly find model fit for different values of alpha:
res.squeezy2 <- squeezy(Y,X, #observed data and response to train model
                        groupset=groupset, #informative co-data groupset
                        Y2=Y2,X2=X2, #test data
                        lambdas = res.squeezy$lambdaMR, #fix lambdas at multi-ridge estimate
                        model="linear", #type of model for the response
                        alpha=0.9) #elastic net mixing parameter


  #Select single-group model or multi-group model based on best mAIC
  res.squeezy <- squeezy(Y,X, #observed data and response to train model
                         groupset=groupset, #informative co-data group set
                         Y2=Y2,X2=X2, #test data
                         fit.ecpc = res.ecpc, #ecpc-fit for initial values
                         model="linear", #type of model for the response
                         alpha=0.5, #elastic net mixing parameter
                         selectAIC = TRUE,resultsAICboth = TRUE)
  
  res.squeezy$modelbestAIC #selected model
  res.squeezy$AICmodels$multigroup$fit$MSEApprox #MSE on test set of multi-group model
  res.squeezy$AICmodels$onegroup$fit$MSEApprox #MSE on test set of single-group model