Q:

Multinomial probit regression with mixed type explanatory variables

I have a data frame called aggregates composed of numerical columns each having a significant amount of zero values. I want to fit a probit model for each column by regressing them on another data frame called exp_vars. exp_vars is composed of factors, ordered factors, integers and numbers. I tried this:
probit_reg = lapply(aggregates, function(y) glm(y ~ subset(exp_vars, select=-c(HH_Net_Income)) + log(exp_vars$HH_Net_Income) + 
log(exp_vars$HH_Net_Income)^2 + log(exp_vars$HH_Net_Income)^3 , data = cbind(y = y, exp_vars)))

But I got this error:
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels 

I know that the error is not related to the formula I set, because even with the most basic formula I get the same error: 
probit_reg=lapply(aggregates, function(y) glm(y ~ exp_vars, data = cbind(y = y, exp_vars)))

I feel like this has something to do with the fact that some of my predictors are factors, but I don't know why it should be a problem and how to fix it. Should I create binary variables for each level of all my factors or is there an easier why to use factor predictors with glm? 
Here is a reproducible example:
    aggregates <- data.frame(agg1 =c(1,0,0,0,0), agg2=c(1,1,1,0,1))
    exp_vars <- data.frame(exp1 =c(21:25), exp2=c(11,12,21,22,23), exp3=c(1,0,0,1,1))
    exp_vars$exp2 = ordered(exp_vars$exp2, levels = c('11','12','21','22','23'))
    exp_vars$exp3 = as.factor(exp_vars$exp3)

     agg1 agg2
#1    1    1
#2    0    1
#3    0    1
#4    0    0
#5    0    1

# exp1 exp2 exp3
#1   21   11    1
#2   22   12    0
#3   23   21    0
#4   24   22    1
#5   25   23    1

probit_reg=lapply(aggregates , function(y) glm(y ~ exp_vars ,family = binomial(link = "probit"), data = cbind(y = y, exp_vars)))

#Error in model.frame.default(formula = y ~ exp_vars , data = cbind(y = y, exp_vars ),  : 
  invalid type (list) for variable 'exp_vars '

A:

The problem is indeed about the formula. The following works
lapply(aggregates , function(y) glm(y ~ . ,family = binomial(link = "probit"), data = cbind(y = y, exp_vars)))

The issue is that when you write y ~ exp_vars, glm looks for a variable named exp_vars and finds a data frame, which it doesn't like. Writing
data = cbind(y = y, exp_vars)

doesn't mean that all the columns of exp_vars can be referred to as exp_vars. In fact,
cbind(y = 1:5, exp_vars)
#   y exp1 exp2 exp3
# 1 1   21   11    1
# 2 2   22   12    0
# 3 3   23   21    0
# 4 4   24   22    1
# 5 5   25   23    1

Meanwhile, writing y ~ . means that regress y on everything else that you find in data, which is what you want.

