Why am I getting such poor results from this example of LASSO-regression? I'm modelling the number of people in a population working, married, being ill, studying, having kids etc. I call one combination of these Boolean variables a "state". 
I'm estimating transition probabilities between these states from one year to the next year. I want to use a logistic regression model to estimate these transition probabilities. Since there are a lot of variables in the data set, I want to try to choose variables automatically and I'm trying LASSO-regression at the moment. 
Below is a made up example where I try to model the probability of working next year. 
The usual glm-function in R works fine and gives estimations of the true parameters that are quite close (although the true influence of the variable married is actually 0). However, when I try LASSO-regression I get results that are far away from the true parameters. Furthermore, I thought that the LASSO-regression would pick up that the variable marriage does not have any influence, but keep the rest of the parameters. Why doesn’t it work?
library(glmnet)
library(boot)

seed <- 1

# Create the combinations of variables
set_working <- 0:1
set_married <- 0:1
set_age <- seq(from=20, to=60, by=1)

data_prob <- expand.grid(set_age, set_working,set_married)
colnames(data_prob) <- c("age","working","married")

# Create true underlying prob-model
intercept_parameter <- -3
age_parameter <- 1/20
working_parameter <- 4
married_parameter <- 0
true_parameters <- c(intercept_parameter,age_parameter,working_parameter,married_parameter)

data_prob$prob_working_next_year <- inv.logit(intercept_parameter + 
                                                data_prob$age * age_parameter + 
                                                data_prob$working * working_parameter +
                                                data_prob$married * married_parameter)

# Generate random observations from the model
set.seed(seed)
data_prob$obs_working_next_year <- rbinom(n = nrow(data_prob), size = 1, prob=data_prob$prob_working_next_year)

# Fit glm-model using all the variables
fit_glm_all_variables <- glm(obs_working_next_year ~ age + working + married  ,
                          family = binomial,
                          data = data_prob)
# The result is reasonably close to the true parameters
coefficients(fit_glm_all_variables)

# Fit glm-model using all the variables and the LASSO-tecnique
x_matrix <- as.matrix(cbind(data_prob$age,data_prob$working,data_prob$married))

fit_glm_lasso <- cv.glmnet(x=x_matrix, 
                           y=data_prob$obs_working_next_year,
                           family = "binomial",
                           alpha=1)

# The result is far away from the true parameters
coefficients(fit_glm_lasso)
# 4 x 1 sparse Matrix of class "dgCMatrix"
# 1
# (Intercept) -0.5932176
# V1           .        
# V2           2.6637460
# V3           .        

# Get the whole sequence of parameters from using different lambda-values
fit_glm_lasso2 <- glmnet(x=x_matrix, 
                         y=data_prob$obs_working_next_year,
                         family = "binomial",
                         alpha=1)
coefficients(fit_glm_lasso2)

# Get the whole sequence of parameters from using different lambda-values
# This time explicately setting the glmnet-function to do the
# standardization
fit_glm_lasso3 <- glmnet(x=x_matrix, 
                         y=data_prob$obs_working_next_year,
                         family = "binomial",
                         alpha=1,
                         standardize = TRUE,
                         intercept = TRUE,
                         standardize.response = FALSE)
coefficients(fit_glm_lasso3)

# The "best" set of parameters chosen by the cv.glmnet-function
# When setting the glmnet-function to do the
# standardization
fit_glm_lasso4 <- cv.glmnet(x=x_matrix, 
                         y=data_prob$obs_working_next_year,
                         family = "binomial",
                         alpha=1,
                         standardize = TRUE,
                         intercept = TRUE,
                         standardize.response = FALSE)
coefficients(fit_glm_lasso4)
# 4 x 1 sparse Matrix of class "dgCMatrix"
# 1
# (Intercept) -0.6273374
# V1           .        
# V2           2.7786131
# V3           .        
# > 

 A: One thing is that your data isn't scaled; this means that age will have an artificially low coefficient compared to married. It is also why age is what drops out first, despite its true influence being higher than married.
A good diagnostic is to plot the path for all lambda:
fit_glm_lasso <- glmnet(x=x_matrix, 
                        y=data_prob$obs_working_next_year,
                        family = "binomial",
                        alpha=1)
plot(fit_glm_lasso)

where the scaling problem is apparent. If you scale your data first, then the lasso will drop out married first, as predicted:
x_scaled = scale(x_matrix)
fit_glm_lasso <- glmnet(x=x_scaled, 
                        y=data_prob$obs_working_next_year,
                        family = "binomial",
                        alpha=1)
plot(fit_glm_lasso)

The coefficients need to be reverse-transformed (not familiar with R enough to do this without a hacky solution), but they are ultimately much closer to your true model, especially if you take the point on the path where married drops out, before the optimum CV value is found.
