```{r}
pacman::p_load(
  tidyverse,       
  glmmLasso)  
```


Data preparation 

```{r}
#  load data

daten = read_rds("freelive.rdata")


daten = daten[complete.cases(daten),]

```

scaling before modelling

```{r}
mets <- daten

mets[, c(3:1149)] <- scale(mets[, c(3:1149)], center = TRUE, scale = TRUE)
```


Define variables 

```{r}
#  specify dependent and group variable

y = mets$WGRye

id = mets$ID

#  specify independent variables

X <- mets[, 3:1149]
```

Define design matrices

```{r}

X = X
Z = matrix(1, length(id), 1)

```


Model computation: Construct data set for glmmLasso 

```{r}
names = c()
for(i in 1:ncol(X)){
  names = c(names, paste0('V', i))
}
colnames(X) = names

mets = as.data.frame(cbind(y, X, id))
mets$id = as.factor(mets$id)

formel = 'y ~ V1'
for(i in 2:ncol(X)){
  formel = paste0(formel, ' + V', i)
}

```

Quick fit  for testing the model:

```{r}
fit = glmmLasso(as.formula(formel), rnd = list(id = ~ 1), lambda=30, data = mets,
                family = gaussian(link="identity"))

```


MAIN CODE for selecting best BIC model
#-----------------

t.las = system.time({
  LAMBDA = seq(1, 30, .1)
  ic = coef = c()
  delta.start = as.matrix(t(rep(0, 1 + ncol(X) + length(unique(id)))))
  delta.start[1, 1] = mean(y)
  Q.start = list()
  Q.start[[1]] = .1
  for(i in 1:length(LAMBDA)){
    print(i)
    mod = tryCatch(glmmLasso(as.formula(formel), rnd = list(id = ~ 1), data = mets, lambda = LAMBDA[i], control = list(start = delta.start[i,], q_start = Q.start[[i]])), error = function(x) Inf)
    if(is.list(mod)){
      ic = c(ic, mod$bic)
      coef = cbind(coef, mod$coefficients)
      delta.start = rbind(delta.start, mod$Deltamatrix[mod$conv.step,])
      Q.start[[i+1]] = as.numeric(mod$Q_long[[mod$conv.step + 1]])
    }else{
      ic = c(ic, Inf)
      coef = cbind(coef, coef[,i-1])
      delta.start = rbind(delta.start, delta.start[i,])
      Q.start[[i+1]] = as.numeric(Q.start[[i]])
    }
  }
  best = which.min(ic)
  
  m.las = glmmLasso(as.formula(formel), rnd = list(id = ~ 1), data = mets, lambda = LAMBDA[i], control = list(start = delta.start[best,], q_start = Q.start[[best]]))
})

#---------------------------------------

cf.las = m.las$coefficients

#----------------------------------------


Results
------------------------------------------------------------------------------------------------
Error after running for a while: 

Error in est.glmmLasso.RE(fix = fix, rnd = rnd, data = data, lambda = lambda, : 
  Fisher matrix not invertible


the error could be found on the original documentation of the package: 

https://rdrr.io/cran/glmmLasso/src/R/glmm_final.smooth_noRE.R


#Trace back results leading to the error: 
--------------------------------------------------------------------------------------------------
4.
stop("Fisher matrix not invertible")
3.
est.glmmLasso.RE(fix = fix, rnd = rnd, data = data, lambda = lambda, 
family = family, final.re = final.re, switch.NR = switch.NR, 
control = control)
2.
est.glmmLasso(fix, rnd, data = data, lambda = lambda, family = family, 
switch.NR = switch.NR, final.re = final.re, control = control)
1.
glmmLasso(as.formula(formel), rnd = list(id = ~1), lambda = 30, 
data = mets, family = gaussian(link = "identity"))