Call:
glm(formula = k ~ a + b + c + d + e + f + g + h + i + j,
    family = binomial(link = "logit"), data = bz)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-2.33397  -1.11978   0.06133   1.12395   2.47743

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)
(Intercept)        -0.442743   0.037655 -11.758  < 2e-16 ***
a                  -0.042182   0.002231 -18.911  < 2e-16 ***
b                   0.514025   0.037674  13.644  < 2e-16 ***
c                  -2.640015   0.166331 -15.872  < 2e-16 ***
d                   1.505434   0.090759  16.587  < 2e-16 ***
e                   1.503102   0.096854  15.519  < 2e-16 ***
f                  -1.262869   0.116334 -10.856  < 2e-16 ***
g                   0.745737   0.179957   4.144 3.41e-05 ***
h                   0.312694   0.021166  14.774  < 2e-16 ***
i                  -0.440660   0.032558 -13.535  < 2e-16 ***
j                   0.773453   0.036602  21.131  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 83178  on 59999  degrees of freedom
Residual deviance: 79092  on 59989  degrees of freedom
AIC: 79114

Number of Fisher Scoring iterations: 4

probabilities <- mod %>% predict(bz, type = "response")
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")

prop.table(table(predicted.classes,bz$k))

predicted.classes         0         1
                0 0.3039000 0.1946333
                1 0.1961000 0.3053667

mean(predicted.classes == bz$millennials_01)
[1] 0.6092667