<p>For lasso+FE, you can first demean both sides of your regression by following the logic given e.g. <a href="https://en.wikipedia.org/wiki/Fixed_effects_model#Classical_Representation" rel="nofollow noreferrer">here</a>, and then run lasso via glmnet. </p>

<p>Lasso+random effects is a bit <a href="https://arxiv.org/abs/1109.4003" rel="nofollow noreferrer">more complicated beast</a> mathematically and it is not supported out of the box with glmnet. There exists a package for doing a mixed-model lasso <a href="https://cran.r-project.org/web/packages/glmmLasso/index.html" rel="nofollow noreferrer">here</a>, but I haven't tried it.</p>
Panel Lasso Estimation in R
<p>I have many continuous independent variables and a dependent dummy variable in my data set about individuals in given years. I want to perform feature selection using Logistic Random Effects Lasso/Logistic Fixed Effects Lasso. However, the default settings of <code>glmnet</code> for my estimation procedure is that I am using cross-sectional data while I want <code>R</code> to see my data as panel data, and it thus models a Logistic Lasso while I want a Logistic Random Effects Lasso/Logistic Fixed Effects Lasso model.</p>

<p>Therefore, in the example code below, I want to let <code>R</code> know that I am using a panel data set and that <code>ID</code> are my individuals/cross-sectional units and <code>year</code> are the years I have observations for each <code>ID</code>. In the code below, all individuals are pooled and I even get coefficients for <code>ID</code> (and <code>year</code>) in this Logistic Lasso estimation. How can I estimate a Logistic Random Effects Lasso/Logistic Fixed Effects Lasso model in <code>R</code>?</p>

<pre><code>df=cbind(c(1,546,2,56,6,73,4234,436,647,567,87,2,5,76,5,456,6756,6,132,78,32),c(2,3546,26,568,76,873,234,36,67,57,887,29,50,736,51,56,676,62,32,782,322),10:30)
year=rep(1:3, times=7)
ID=rep(1:7, each=3)
x=as.matrix(cbind(ID,year,df))
y1=as.data.frame(rep(c(0,1), each = 18))[1:21,]
y=as.matrix(y1)

fit=glmnet(x,y,alpha=1,family="binomial")
lambdamin=min(fit$lambda)
predict.glmnet(fit,s=lambdamin,newx=x,type="coefficients")
</code></pre>

<pre><code>                        1
(Intercept) -8.309211e+01
ID           1.281220e+01
year         .           
            -2.339904e-04
             .           
             .           
</code></pre>