sAIC | R Documentation |
This function computes the Akaike information criterion for generalized linear models estimated by the lasso.
sAIC(x, y=NULL, beta, family=c("binomial","poisson","ggm"))
x |
A data matrix. |
y |
A response vector. If you select family="ggm", you should omit this argument. |
beta |
An estimated coefficient vector including the intercept. If you select family="ggm", you should use an estimated precision matrix. |
family |
Response type (binomial, Poisson or Gaussian graphical model). |
AIC |
The value of AIC. |
Shuichi Kawano
skawano@math.kyushu-u.ac.jp
Ninomiya, Y. and Kawano, S. (2016). AIC for the Lasso in generalized linear models. Electronic Journal of Statistics, 10, 2537–2560. doi: 10.1214/16-EJS1179
library(MASS) library(glmnet) library(glasso) ### logistic model set.seed(3) n <- 100; np <- 10; beta <- c(rep(0.5,3), rep(0,np-3)) Sigma <- diag( rep(1,np) ) for(i in 1:np) for(j in 1:np) Sigma[i,j] <- 0.5^(abs(i-j)) x <- mvrnorm(n, rep(0, np), Sigma) y <- rbinom(n,1,1-1/(1+exp(x%*%beta))) glmnet.object <- glmnet(x,y,family="binomial",alpha=1) coef.glmnet <- coef(glmnet.object) ### coefficients coef.glmnet[ ,10] ### AIC sAIC(x=x, y=y, beta=coef.glmnet[ ,10], family="binomial") ### Poisson model set.seed(1) n <- 100; np <- 10; beta <- c(rep(0.5,3), rep(0,np-3)) Sigma <- diag( rep(1,np) ) for(i in 1:np) for(j in 1:np) Sigma[i,j] <- 0.5^(abs(i-j)) x <- mvrnorm(n, rep(0, np), Sigma) y <- rpois(n,exp(x%*%beta)) glmnet.object <- glmnet(x,y,family="poisson",alpha=1) coef.glmnet <- coef(glmnet.object) ### coefficients coef.glmnet[ ,20] ### AIC sAIC(x=x, y=y, beta=coef.glmnet[ ,20], family="poisson") ### Gaussian graphical model set.seed(1) n <- 100; np <- 10; lambda_list <- 1:100/50 invSigma <- diag( rep(0,np) ) for(i in 1:np) { for(j in 1:np) { if( i == j ) invSigma[i ,j] <- 1 if( i == (j-1) || (i-1) == j ) invSigma[i ,j] <- 0.5 } } Sigma <- solve(invSigma) x <- scale(mvrnorm(n, rep(0, np), Sigma)) glasso.object <- glassopath(var(x), rholist=lambda_list, trace=0) ### AIC sAIC(x=x, beta=glasso.object$wi[,,10], family="ggm")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.