estimate: Penalized lvm model

In bozenne/lava.penalty: Penalized Latent Variable Models

Description

Estimate a penalized lvm model

Usage

## S3 method for class 'plvm'
estimate(x, data, estimator = "gaussian1", lambda1, lambda2,
  lambdaG, lambdaN, adaptive = FALSE, regularizationPath = FALSE,
  fit = lava.options()$calcLambda$fit, control = list(),
  control.proxGrad = list(), control.EPSODE = list(),
  equivariance = FALSE, ...)

## S3 method for class 'pmultigroup'
estimate(x, data2, estimator = "gaussian1", lambda1,
  lambda2, lambdaG, lambdaN, adaptive = FALSE, regularizationPath = FALSE,
  fit = lava.options()$calcLambda$fit, control = list(),
  control.proxGrad = list(), control.EPSODE = list(),
  equivariance = FALSE, ...)

Arguments

x	a penalized lvm model
data	a data.frame containing the data
estimator	the method used to compute the likelihood and its first two derivatives.
lambda1	lasso penalization parameter
lambda2	ridge penalization parameter
lambdaG	group lasso penalization parameter
lambdaN	nuclear norm penalization parameter
adaptive	should the coefficient of the adaptive lasso be returned instead of the coefficient of the lasso
regularizationPath	should the regularization path be computed. If so the argument lambda1 is ignored but not lambda2.
fit	the criterion used to retain the optimal lambda1 value when the regularization path is used.
control	control/optimization parameters
control.proxGrad	control values for the proximal gradient algorithm (see proxGrad)
control.EPSODE	control values for the EPSODE algorithm (see EPSODE)
equivariance	should the regularization parameters be multiplied by the first variance parameter and the first variance parameter fixed to 0.
...	additional arguments to be passed to lava:::estimate.lvm

Details

Available estimators are:

• "penalized" the hessian is computed as E[t(S) S]

• "numDeriveSimple" the hessian is computed using numerical derivatives

• "numDeriveRichardson" the hessian is computed using numerical derivatives (Richardson method)

• "explicit" the hessian is computed using an explicit formula

Value

a plvmfit object

References

Zhou 2014 - A generic Path Algorithm for Regularized Statistical Estimation
Park 2007 - L1-regularization path algorithm for GLM

Examples

# library(butils.base)
# package.source("lavaPenalty", RorderDescription = FALSE)
library(lavaPenalty)
library(ggplot2)

options.proxGrad <- lava.options()$proxGrad
options.proxGrad$export.iter <- TRUE
lava.options(proxGrad = options.proxGrad)

####> linear regression ####
set.seed(10)
n <- 300
formula.lvm <- as.formula(paste0("Y~",paste(paste0("X",1:12), collapse = "+")))
mSim <- lvm(formula.lvm)
df.data <- sim(mSim,n)

lvm.model <- lvm(formula.lvm)
plvm.model <- penalize(lvm.model)
penalty(plvm.model)
plvm.model$penalty

#### lasso penalty

## regularization path
path1F <- estimate(plvm.model,  data = df.data, regularizationPath = TRUE)
plot(path1F, type = "path")
plot(path1F, type = "path", coefficient = "npenalized")

## proxGrad
pfit <- estimate(plvm.model,  data = df.data,
                 lambda1 = getPath(path1B, names = "lambda1.abs")[6,1],
                 control = list(trace = 2, constrain = TRUE), fixSigma = TRUE)
pfit

#### ridge penalty
pfit <- estimate(plvm.model,  data = df.data, lambda2 = 10, control = list(constrain = TRUE))
pfit

#### elastic net penalty
pfit <- estimate(plvm.model,  data = df.data, lambda1 = 100, lambda2 = 10, control = list(constrain = TRUE))
pfit

#### group lasso penalty
plvm.model <- penalize(lvm.model, value = paste0("Y~X",1:4), group = 1)
plvm.model <- penalize(plvm.model, value = paste0("Y~X",5:9), group = 2)

pfit <- estimate(plvm.model,  data = df.data, lambda1 = 80, control = list(constrain = TRUE), fixSigma = TRUE)
pfit

#### algorithm ISTA
lambda1 <- 10

pfit_ISTA <- estimate(plvm.model,  data = df.data, method.proxGrad = "ISTA",
                      lambda1 = lambda1, iter.max = 100,
                      control = list(constrain = TRUE), fixSigma = TRUE)

pfit_FISTA_Vand <- estimate(plvm.model,  data = df.data, method.proxGrad = "FISTA_Beck",
                            lambda1 = lambda1, iter.max = 100,
                            control = list(constrain = TRUE), fixSigma = TRUE)

pfit_FISTA_Beck <- estimate(plvm.model,  data = df.data, method.proxGrad = "FISTA_Vand",
                            lambda1 = lambda1, iter.max = 100,
                            control = list(constrain = TRUE), fixSigma = TRUE)

pfit_mFISTA_Vand <- estimate(plvm.model,  data = df.data, method.proxGrad = "mFISTA_Vand",
                             lambda1 = lambda1, iter.max = 100,
                             control = list(constrain = TRUE), fixSigma = TRUE)


df.cv <- rbind(data.frame(pfit_ISTA$opt$details.cv, method = "ISTA"),
               data.frame(pfit_FISTA_Vand$opt$details.cv, method = "FISTA_Beck"),
               data.frame(pfit_FISTA_Beck$opt$details.cv, method = "FISTA_Vand"),
               data.frame(pfit_mFISTA_Vand$opt$details.cv, method = "mFISTA_Vand"))
gg <- ggplot(df.cv, aes(x = iteration, y = obj, group = method, color = method)) 
gg <- gg + geom_line() + geom_point()
gg

gg + coord_cartesian(ylim = c(1490,1510), xlim = c(0,100))
gg + coord_cartesian(ylim = c(1500,1510), xlim = c(0,100))

#### nuclear norm penalty
n.obs <- 100
res <- simForm(n.obs, xmax = 25, ymax = 25, radius = 5)
coords <- res$coords
n.coord <- nrow(coords)
betaI <- as.vector(res$X)
X <- matrix(rnorm(n.obs*n.coord), nrow = n.obs, ncol = n.coord)
Xnames <- paste0("X",1:n.coord)

n.confounder <- 5
gamma <- rep(1, n.confounder)
Z <- matrix(rnorm(n.obs*n.confounder), nrow = n.obs, ncol = n.confounder)
Znames <- paste0("Z",1:n.confounder)

Y <- Z %*% gamma + X %*% betaI + rnorm(n.obs)
formula.lvm <- as.formula( paste0("Y~", paste0("X",1:n.coord, collapse = "+") ) )

dt.data <- data.table(Y=Y,data.frame(Z),data.frame(X))
names(dt.data) <- c("Y",Znames,Xnames)
dt.data[, (names(dt.data)) := lapply(.SD,scale), .SDcols = names(dt.data)]

##
lvm.image <- lvm(as.formula(paste0("Y~",paste(Znames, collapse = "+"))))
plvm.image <- lvm.image
penalizeNuclear(plvm.image, coords = coords) <- as.formula(paste0("Y~",paste0(Xnames,collapse = "+")))

start <- setNames(rep(0, length = length(coef(plvm.image))), coef(plvm.image))
start["Y,Y"] <- 1

for(lambda in c(1e0,1e1,1e2,2e2)){
  elvm.Path <- estimate(plvm.image,  data = dt.data, lambdaN = lambda, method.proxGrad = "FISTA_Vand",
                        control = list(trace = 2,iter.max = 10))
  B.LS <- matrix(attr(elvm.Path$opt$message,"par")[paste0("Y~",Xnames)],
                 nrow = NROW(res$X), ncol = NCOL(res$X), byrow = TRUE)
  fields:::image.plot(B.LS)
  
  plot(elvm.Path$opt$details.cv[,"obj"])
  
}


####> Latent variable model ####

set.seed(10)
n <- 300
mSim <- lvm(list(Y1 ~ eta + X1, Y2 ~ eta, Y3 ~ eta + X2))
latent(mSim) <- ~eta
regression(mSim) <- eta~Z1
covariance(mSim) <- Y1~Y2

df.data <- sim(mSim,n)


#### partial lasso penalty
m <- mSim
regression(m) <- Y2~X1+X2
pm <- penalize(m, c("Y2~X1","Y2~X2","Y3~X2"))

## regularization path
path.lvm <- estimate(pm, df.data, regularizationPath = TRUE, stopParam = 1, increasing = FALSE)

pathLVM <- calcLambda(model = pm, data.fit = df.data, seq_lambda1 = seq(0,300, length.out = 5))
plot(pathLVM)
plot(pathLVM, type = "path")
# lava:::estimate.lvm(pm, scale(df.data))

#### algorithm ISTA
lambda1 <- 10

e_ISTA <- estimate(pm,  data = df.data, method.proxGrad = "ISTA",
                      lambda1 = lambda1,
                      control = list(constrain = TRUE))

e_FISTA_Vand <- estimate(pm,  data = df.data, method.proxGrad = "FISTA_Beck",
                            lambda1 = lambda1,
                            control = list(constrain = TRUE))

e_FISTA_Beck <- estimate(pm,  data = df.data, method.proxGrad = "FISTA_Vand",
                            lambda1 = lambda1,
                            control = list(constrain = TRUE))

e_mFISTA_Vand <- estimate(pm,  data = df.data, method.proxGrad = "mFISTA_Vand",
                             lambda1 = lambda1,
                             control = list(constrain = TRUE))


df.cv <- rbind(data.frame(e_ISTA$opt$details.cv, method = "ISTA"),
               data.frame(e_FISTA_Vand$opt$details.cv, method = "FISTA_Beck"),
               data.frame(e_FISTA_Beck$opt$details.cv, method = "FISTA_Vand"),
               data.frame(e_mFISTA_Vand$opt$details.cv, method = "mFISTA_Vand"))
gg <- ggplot(df.cv, aes(x = iteration, y = obj, group = method, color = method)) 
gg <- gg + geom_line() + geom_point() +
gg

bozenne/lava.penalty documentation built on May 13, 2019, 1:41 a.m.