Seguindo a opinião do Daniel rodei um modelo usando Poisson, porém o mesmo não foi adequado (Resíduos e Deviance), não sei se os número elevado de 0 em algumas co-variável possam influenciar nos resultados. Então rodei um modelo Binomial Negativo que teve melhores resultados, como nunca trabalhei com esses modelos estou com algumas dúvidas

Como escrever colocar os resultados do modelo Binomial Negativo em uma função matemática?

Qual o comando que posso usar para predição?

Call:
glm.nb(formula = esqui ~ alfabetizado + Fossa.rudimentar...não.discriminada +
    Coletado.por.caçamba.de.serviço.de.limpeza + Poço.ou.nascente...canaliz.em.pelo.menos.um.cômodo +
    Poço.ou.nascente...canaliz.só.na.propried.terreno + Outra.forma...não.canalizada +
    offset(log(pop)), data = dados, init.theta = 0.5883720786,
    link = log)

Deviance Residuals:
Min 1Q Median 3Q Max
-2.3668 -1.2647 -0.4984 0.1509 2.2684

Coefficients:
                                                   Estimate Std. Error z value Pr(>|z|)
(Intercept)                                         -0.9340     0.7251 -1.288 0.197687
alfabetizado                                        -6.2288     0.9966 -6.250 4.11e-10 ***
Fossa.rudimentar...não.discriminada                  0.6584     0.3802   1.732 0.083276 .
Coletado.por.caçamba.de.serviço.de.limpeza           1.5681     0.7560   2.074 0.038062 *
Poço.ou.nascente...canaliz.em.pelo.menos.um.cômodo   2.2748     0.6546   3.475 0.000511 ***
Poço.ou.nascente...canaliz.só.na.propried.terreno    4.8851     1.7172   2.845 0.004445 **
Outra.forma...não.canalizada                        -1.8832     0.7683 -2.451 0.014239 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.5884) family taken to be 1)

Null deviance: 396.86 on 268 degrees of freedom
Residual deviance: 329.11 on 262 degrees of freedom
AIC: 3379.9

Number of Fisher Scoring iterations: 1

Theta: 0.5884
Std. Err.: 0.0432

2 x log-likelihood: -3363.9160

--