Introduction to gglasso

Install the package

# on CRAN
install.packages("gglasso")

# dev version on GitHub
pacman::p_load_gh('emeryyi/gglasso')

Least squares regression

library(gglasso)

# load bardet data set
data(bardet)

group1 <- rep(1:20, each = 5)

fit_ls <- gglasso(x = bardet$x, y = bardet$y, group = group1, loss = "ls")

plot(fit_ls)

coef(fit_ls)[1:5,90:100]

##                      s89          s90          s91          s92          s93
## (Intercept)  8.099354325  8.098922472  8.098531366  8.098175719  8.097849146
## V1          -0.119580203 -0.120877799 -0.122079683 -0.123223779 -0.124310183
## V2          -0.113742329 -0.114834411 -0.115853997 -0.116837630 -0.117782854
## V3          -0.002584792 -0.003487571 -0.004328519 -0.005134215 -0.005904892
## V4          -0.084771705 -0.088304073 -0.091674509 -0.094978960 -0.098212775
##                      s94          s95          s96          s97          s98
## (Intercept)  8.097574095  8.097295166  8.097058895  8.096833259  8.096637676
## V1          -0.125274109 -0.126284595 -0.127173301 -0.128011016 -0.128738414
## V2          -0.118630121 -0.119526679 -0.120326451 -0.121086672 -0.121754134
## V3          -0.006593702 -0.007323107 -0.007970047 -0.008583011 -0.009116809
## V4          -0.101161988 -0.104349829 -0.107241330 -0.110045942 -0.112543755
##                      s99
## (Intercept)  8.096455264
## V1          -0.129453437
## V2          -0.122415680
## V3          -0.009645386
## V4          -0.115058449

Cross-Validation

cvfit_ls <- cv.gglasso(x = bardet$x, y = bardet$y, group = group1, loss = "ls")

plot(cvfit_ls)

coef(cvfit_ls, s = "lambda.min")

##                         1
## (Intercept)  8.249757e+00
## V1           2.542969e-03
## V2          -4.251746e-02
## V3           2.619704e-02
## V4           7.010593e-03
## V5          -7.853723e-02
## V6           2.056230e-04
## V7           3.494186e-04
## V8           1.861031e-04
## V9           8.471457e-06
## V10         -7.538880e-04
## V11          2.472638e-03
## V12         -1.515146e-03
## V13          5.038118e-05
## V14         -4.159035e-06
## V15         -5.268445e-03
## V16          1.207185e-02
## V17          3.757153e-02
## V18         -1.813877e-02
## V19         -5.866357e-03
## V20         -7.578878e-02
## V21          6.912649e-02
## V22          7.385719e-02
## V23         -5.115425e-02
## V24         -1.021370e-02
## V25         -2.176734e-01
## V26          7.889141e-02
## V27          4.079558e-02
## V28         -1.040149e-02
## V29          2.147319e-03
## V30         -1.588000e-01
## V31          0.000000e+00
## V32          0.000000e+00
## V33          0.000000e+00
## V34          0.000000e+00
## V35          0.000000e+00
## V36          3.093717e-02
## V37         -2.279141e-02
## V38          8.466720e-03
## V39          1.075600e-02
## V40         -6.194654e-02
## V41          0.000000e+00
## V42          0.000000e+00
## V43          0.000000e+00
## V44          0.000000e+00
## V45          0.000000e+00
## V46         -2.320007e-02
## V47          1.218448e-02
## V48          3.821948e-02
## V49          8.642710e-03
## V50          5.080079e-03
## V51         -1.455627e-02
## V52         -1.056218e-02
## V53          8.577112e-02
## V54          1.333652e-01
## V55          3.847231e-02
## V56          0.000000e+00
## V57          0.000000e+00
## V58          0.000000e+00
## V59          0.000000e+00
## V60          0.000000e+00
## V61         -3.006023e-02
## V62          7.142155e-03
## V63          5.018859e-03
## V64          4.961125e-02
## V65          9.332342e-02
## V66         -1.280058e-02
## V67          7.584048e-04
## V68          2.086525e-02
## V69          1.758820e-02
## V70          1.422012e-03
## V71         -7.634636e-03
## V72         -5.527457e-03
## V73          2.479708e-02
## V74          4.632734e-06
## V75          1.429469e-02
## V76         -9.886164e-03
## V77          6.116407e-03
## V78          2.597462e-02
## V79          1.266998e-03
## V80          1.012356e-02
## V81          4.084980e-03
## V82         -6.410604e-03
## V83          2.612768e-03
## V84          2.789949e-03
## V85         -1.316906e-02
## V86          1.171404e-02
## V87         -3.658167e-02
## V88          1.193755e-02
## V89          1.906166e-02
## V90         -5.838967e-02
## V91          0.000000e+00
## V92          0.000000e+00
## V93          0.000000e+00
## V94          0.000000e+00
## V95          0.000000e+00
## V96          0.000000e+00
## V97          0.000000e+00
## V98          0.000000e+00
## V99          0.000000e+00
## V100         0.000000e+00

Weight Least squares regression

We can also perform weighted least-squares regression by specifying loss='wls', and providing a \(n \times n\) weight matrix in the weights argument, where \(n\) is the number of observations. Note that cross-validation is NOT IMPLEMENTED for loss='wls'.

# generate weight matrix
times <- seq_along(bardet$y)
rho <- 0.5
sigma <- 1
H <- abs(outer(times, times, "-"))
V <- sigma * rho^H
p <- nrow(V)
V[cbind(1:p, 1:p)] <- V[cbind(1:p, 1:p)] * sigma

# reduce eps to speed up convergence for vignette build
fit_wls <- gglasso(x = bardet$x, y = bardet$y, group = group1, loss = "wls", 
                   weight = V, eps = 1e-4)

plot(fit_wls)

coef(fit_wls)[1:5,90:100]

##                     s89         s90         s91         s92         s93
## (Intercept)  8.09429262  8.09340481  8.09254573  8.09170743  8.09089247
## V1          -0.13922372 -0.14077803 -0.14222609 -0.14359110 -0.14487482
## V2          -0.15966042 -0.16117772 -0.16261019 -0.16397683 -0.16527730
## V3           0.03917529  0.03880296  0.03847035  0.03816594  0.03788642
## V4          -0.16548208 -0.17057112 -0.17546237 -0.18021267 -0.18481370
##                     s94         s95         s96         s97         s98
## (Intercept)  8.09011527  8.08935394  8.08862146  8.08793054  8.08727098
## V1          -0.14606257 -0.14719352 -0.14824987 -0.14921624 -0.15011520
## V2          -0.16649386 -0.16766459 -0.16877053 -0.16979410 -0.17075592
## V3           0.03763241  0.03739356  0.03717453  0.03697953  0.03680099
## V4          -0.18919074 -0.19347289 -0.19758499 -0.20145156 -0.20513769
##                     s99
## (Intercept)  8.08664325
## V1          -0.15094837
## V2          -0.17165664
## V3           0.03663899
## V4          -0.20863833