Beni
Beni

Reputation: 11

Use of glm with ri2

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

Answers (0)

Related Questions