Reputation: 11
I am conducting a randomisation inference procedure with ri2 for a binary model. I use glm and specify as a family binomial (link = 'logit'). When I use the ri2 on this model apparently it is not considering this family. I get different coefficients which seem to be similar to the output when not specifying any family; the default family seems to be gaussian which equates to using lm if I am understanding this correctly.
I was wondering then whether in this case this randomisation inference test can be considered valid for the model that I specify? At the end I would like to use the ri2 as a check on the significance level of the treatment arm variable with the glm by comparing the p-value of the coefficients of my treatment arms.
Is maybe there is a way to still have the family correctly included in ri2? I am not sure why the function is not recognised in this package.
Here is my code to illustrate the problem. I compare the results of the glm outputs with and without specified family, and the result of ri2.
library(dplyr)
library(lmtest)
library(ri2)
set.seed(123)
n <- 100
y <- rbinom(100, size = 1, prob = 0.5)
ngroup = 4
nrep = 15
treat <- rep(c("control", "treat1", "treat2", "treat3"), each = nrep)
strata <- sample(1:3, 100, replace =TRUE)
cluster <- sample(1:5, 100, replace =TRUE)
test_1 <- as.data.frame(cbind(treat, strata, cluster, y))
test_1 <- test_1 %>%
mutate(y = as.numeric(y),
cluster = as.numeric(cluster))
test_1 <- fastDummies::dummy_cols(test_1, select_columns = "treat")
test_1 <- fastDummies::dummy_cols(test_1, select_columns = "strata")
### glm with family = binomial(link = "logit")
m1 <- glm(y~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3, data = test_1, family = binomial(link = "logit"))
### glm with a default family
m2 <- glm(y~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3, data = test_1)
### results with ri2
declaration <- declare_ra(N = 100, clusters = test_1$cluster)
ri2 <-
conduct_ri(
formula = glm(y~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3, data = test_1, family = binomial(link = "logit")), # here also only the saved model to-be-tested, m1, can be specified instead of the formula.
assignment = "treat_treat1",
declaration = declaration,
data = test_1,
sims = 1000
)
ri2
coeftest(m1)
coeftest(m2)
Upvotes: 0
Views: 47