American-Institutes-for-Research / WeMix

WeMix public repository
GNU General Public License v2.0
10 stars 2 forks source link

Predict with WeMix #12

Closed Gabatsarl closed 10 months ago

Gabatsarl commented 1 year ago

Hello, I would like to use the "WeMix" package. I would like to know if there is a prediction function in the R software? How can I get the residuals out of the model?

my model is

model2019_mix_lect1 <- mix(LECTURE1 ~ qe21+qe22+qe25+qe28+(1|ID_ECOLE), 
                        data = data2019,weights=c("rwgt0", "W_SCH_ADJ2"),family = "binomial")

summary(model2019_mix_lect1)
# Call:
# mix(formula = LECTURE1 ~ qe21 + qe22 + qe25 + qe28 + (1 | ID_ECOLE), 
#     data = data2019, weights = c("rwgt0", "W_SCH_ADJ2"), family = "binomial")
# 
# Variance terms:
#  Level    Group        Name Variance Std. Error Std.Dev.
#      2 ID_ECOLE (Intercept)    3.447     0.4522    1.857
# Groups:
#  Level    Group n size mean wgt sum wgt
#      2 ID_ECOLE    216    42.54    9188
#      1      Obs   3823    50.22  192006
# 
# Fixed Effects:
#             Estimate Std. Error t value
# (Intercept)  5.14748    0.50548  10.183
# qe21        -0.26734    0.03797  -7.040
# qe22        -0.25841    0.10441  -2.475
# qe25        -0.13653    0.03599  -3.794
# qe28        -0.07406    0.02233  -3.316
# 
# lnl= -87358.61

I know that with package "glmer" it's easy with the "predict" function and the "resid" function for the residuals. But how do you do it with the "WeMix" package? Thanks, Talagbe [updated by PB to format code]

pdbailey0 commented 1 year ago

@Gabatsarl just to pin down exactly what you are asking for, are you looking for output like this?

# example from ?glmer examples
require(lme4)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)
predict(gm1)
#          1          2          3          4          5          6          7 
# -0.8087134 -1.8006384 -1.9369296 -2.3884588 -1.6974362 -2.6893612 -2.8256524 
resid(gm1)
#           1           2           3           4           5           6           7           8 
# -1.43770782  0.98847195  2.35668826 -0.93702271 -0.24317316 -0.14282938 -0.17032930  0.95457653 
Gabatsarl commented 1 year ago

Yes exactly @pdbailey0, but I want it with the model using WeMix.

I have two questions

Question 1: if

model<- mix(READ1 ~ qe21+qe22+qe25+qe28+(1|ID_ECOLE), 
                        data = data2019,weights=c("rwgt0", "W_SCH_ADJ2"),family = "binomial")

can we do predict(model)?

Question 2: for a new dataset newdata, can we also do predit(model,newdata=newdata)? Thanks [edited by PB for format and type in name "WeMix"]

pdbailey0 commented 1 year ago

There is no such function in WeMix. You could request it.

Gabatsarl commented 1 year ago

Thanks. Conclusion: You can't predict with a model fitted with mix.

pdbailey0 commented 1 year ago

Here is some demo code

require(lme4)
require(WeMix)
cbpp2 <- cbpp1 <- cbpp[,c("herd", "period")]
cbpp1$n <- cbpp$size - cbpp$incidence
cbpp1$level <- 0

cbpp2$n <- cbpp$incidence
cbpp2$level <- 1
cbpp_ <- rbind(cbpp1, cbpp2)
repn <- function(x) {
  x[rep(1, x$n),]
}
cbpp_rep <- do.call(rbind, lapply(split(cbpp_, factor(1:nrow(cbpp_))), repn))
cbpp_rep$one <- 1

gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)

gm2 <- glmer(level ~ period + (1 | herd),
               data = cbpp_rep, family = binomial)

we1 <- mix(level ~ period + (1|herd), data=cbpp_rep, family=binomial(), weights=c("one", "one"))
predict(we1) #  does not work
pdbailey0 commented 10 months ago

Thank you for pointing this out @ Gabatsarl. Our returns lacked some important information for using the results and now have that information. In particular, predict does now work. It also works with new data.