model <- glm(income ~ education + age + married + race + left + quartile_health, family=binomial, data=data)

logitmfx(income ~ education + age + married + race + left + quartile_health, data = data,  atmean = TRUE)

coef(r_3)[1] + coef(r_3)[2]*mean_income + coef(r_3)[3]*mean_age + coef(r_3)*mean_married[4] + coef(r_3)[5]*mean_race + coef(r_3)[6] + coef(r_3)[7]*meanhealth