issues with random effect + question output get_pars()

[ValentinJourne]

Heya, thanks for your nice package, I have tried to fit a model with a random effect, but I got an eror. I also tried to use the example you provided (which is weid, because in both case, with my data, or your example, all row got a value).

set.seed(123)
y <- matrix(runif(99,1,4), ncol=3)
design <- data.frame("fac" = sample(letters[1:5], size=nrow(y), replace=TRUE))
design$fac <- as.factor(design$fac)
fit<-fit_zoid(formula=~(1|fac),design_matrix=design,data_matrix=y,chains=1,iter=100)
Error in na.fail.default(list(`1 | fac` = c(NA, NA, NA, NA, NA, NA, NA,  : 
  missing values in object
In addition: Warning message:
In Ops.factor(1, fac) : ‘|’ not meaningful for factors

---

I also wanted to know, since I am here, once you are using get_pars(), in the output $betas, how can I recognize my variable? For example, if I have Y = cbind(var1, var2, var3), when I am looking in the $betas from get_pars function if m = 1, is it my var1 or var2 or var 3? how could we determine the different m? (e.g. which variable is used a reference)

Thank you!

[ValentinJourne]

Ok, so I do not know if it is the best, but from your function fit_zoid, I commented on those lines (because parse_re_formula will already generate the matrix, and extract random effect factor? L16)

#model_frame <- model.frame(formula, design_matrix)
#model_matrix <- model.matrix(formula, model_frame)
res <- parse_re_formula(formula, design_matrix)

[ericward-noaa]

Great questions and thanks for the help in catching issues. On the parameters:

The betas represent the mean estimates of the parameters in logit space, before being back-transformed. Maybe a clearer example would be to use a case of 3 observations (rows) and 4 groups:

y <- matrix(c(3.77, 6.63, 2.60, 0.9, 1.44, 0.66, 2.10, 3.57, 1.33, 3.4, 1.4, 5),
  nrow = 3, byrow = TRUE
)
fit <- fit_zoid(data_matrix = y)
p_hat <- get_pars(fit)

Now, p_hat$betas contains the estimated intercepts in logit space (the last is set to 0 for identifiability)

And then the back-transformed estimates of proportions are given by p_hat$p

Note that the intercepts are not varying across observations -- so they're the same for all groups. More complicated models with fixed / random effects could make them vary

[ericward-noaa]

For the question about the error and random effects, I'm not sure there's anything wrong here? I was able to run the example just fine (I realize there was the warning message, and just have pushed a commit to suppress that message)

set.seed(123)
y <- matrix(runif(99,1,4), ncol=3)
design <- data.frame("fac" = sample(letters[1:5], size=nrow(y), replace=TRUE))
design$fac <- as.factor(design$fac)
fit<-fit_zoid(formula=~(1|fac),design_matrix=design,data_matrix=y,chains=1,iter=100)

What's the results of sessionInfo()? Here's mine:

[ValentinJourne]

Thanks for the fast answer! :) but in your example, how is the variable defined at 0 (the group you use by default)? In your first screenshot we can see that m = 4, the coefficient is at 0. Does this mean that you are using the 4th column from Y as a reference? Or is it the first column? If I am asking, it's because by default in brms::, when you are fitting a dirichlet model, it will use the first column as "reference category"

[ericward-noaa]

Yes -- the last column is used as the baseline / reference. The choice here is arbitrary, and the back transformed proportions will be identical, no matter which column is set to 0

[ValentinJourne]

OK, I found the issue,...! I used in my code

options(na.action = "na.fail")

And because this function was loaded at the beginning of the code, I did not pay enough attention to it. Thanks a lot for your answers and your help!

[ericward-noaa]

Ok great, closing this issue -- thanks!