Where to start: own simulations on logistic regression. Hosmer-Lemeshow, Farrington, Pearson Chi-Square [closed]

Question

I'd like to run my own simulations on the logistic regression model. I want to test on different grouping strategies, spare data and strategies for combined small groups.

I read a lot of papers by Hosmer and Lemeshow refering to logistic regression. Even further papers by Kuß, Allison, Paul, McGullgh etc. Now I'd like to check some results by my own simulations. Problem is, I don't have much experiences with simulations. Can I run my simulations on R?`Is there a better application (SPSS, MATLAB, Python, ... ?)

Is there a guide for such simulations?

Answer 1

A couple of things you need to know before running your own simulations in R:

• random variable generation: see the R probability distributions page https://cran.r-project.org/web/views/Distributions.html
• one or more data generation processes for the model you're trying to test, regardless of software
• process repetition, the simplest being replicate()

So for example, a data generation process for the logistic regression model is something like: $$Y = \mathrm{Bern}\big((1 + e^{-X\beta})^{-1}\big)$$ where $Y$ is an $n$ by 1 vector of 1s and 0s, $X$ is a matrix of $n$ rows and $p$ predictors (including the intercept, so the first column of $X$ is $1_n$), $\beta$ represents $p$ regression coefficients.

So we can build this example from the ground up with a sample size of 200. Let's assume we have two predictor variables, $x_c \sim N(0, 1.5)$ and $x_b \sim \mathrm{Bern}(0.5)$. This means $x_c$ is normally distributed with mean zero and standard deviation 1.5, and $x_b$ is Bernoulli with 50% probability. In R, this would be:

n <- 200
xb <- rbinom(n = n, size = 1, prob = .5) # see ?rbinom to learn details
xc <- rnorm(n = n, mean = 0, sd = 1.5)

Now let's assume $$X\beta=1+0.5 \times x_c - 0.75 \times x_b$$

That is:

Xbeta <- 1 + 0.5 * xc - 0.75 * xb

Finally, to transform this all to 0s and 1s given the first equation above:

y <- rbinom(length(Xbeta), 1, plogis(Xbeta)) # plogis is equal to 1/(1+e^(-Xbeta))

The simplest exercise to do would be to run the model and save the coefficients to check for unbiased estimation:

coef(glm(y ~ xb + xc, binomial))

And we can wrap it all up using the replicate() function. We'll replicate this 2500 times:

set.seed(12345) # do this to obtain reproducible results
n <- 200
results <- replicate(2500, {
  xb <- rbinom(n, 1, .5)
  xc <- rnorm(n, 0, 1.5)
  y <- rbinom(n, 1, plogis(1 + 0.5 * xc - 0.75 * xb))
  coef(glm(y ~ xb + xc, binomial))
})

This will save the three coefficients 2500 times. We'll get 3 rows and 2500 columns. It's nicer to transpose the matrix, so:

results <- t(results)
colMeans(results) # obtain the mean estimated coefficients
(Intercept)          xb          xc
  1.0232931  -0.7790580   0.5151716

You can play around with saving more than the coefficients. Depending on your proficiency with R, you can develop better ways to manage the process. Once you find this barebones approach becoming unwieldy, look up the SimDesign package. In this paper (https://amstat.tandfonline.com/doi/full/10.1080/10691898.2016.1246953), the authors demonstrate its use, so it should not be too difficult to get started.

Personally, I use replicate() to check things very quickly. I use SimDesign when I want additional features that make managing the process easy.