I am trying to simulate possible sample sizes for a proportional odds model using the posamsize function available in the Hmisc package. The model will look at a Y (dependant variable) that has 11 possible levels. After a round of going around I finally came across a simple way to generate a probability distribution for these 11 levels (probability of the number of observations for each level) using the Dirichlet distribution (stackexchange answer). The code I am using is as follows
library(Hmisc)
library(tidyverse)
library(gtools)
df <- as.data.frame(rdirichlet(2000,alpha=c(1:11)))
df <-df %>%
rowwise(.) %>%
mutate(samplesize=posamsize(c_across(V1:V11),odds.ratio=2.0,power=0.8,alpha=0.05)$n)
This generates a data frame with 11 columns and 2000 rows where each column corresponds to the levels of the Y. The row sums are always 1. The sample size using this method seems to range between 198 - 205. This implies that the method is generating probabilities at each level which are quite narrow (for example in V11 the range is between 0.04 - 0.34 and V1 is 0 - 0.13)
I have previously tried a naive method using expand.grid to accommodate all possible combinations of probabilities ranging between 0 - 0.7 at 0.1 increments but ran out of memory.
The method does demonstrate what the BBR notes state that the frequency of observations in each category of Y does not influence that sample size that much - I have noted the largest sample size occurs when there are extreme distributions - e.g. one category gets 95% and another one gets 5% and rest are 0.
I would like to know if the method of simulation I am adopting is correct and if a better method is available that will show sample size for the extreme probability distributions of Y
Edit I modified the simulation code as below:
set.seed(100)
df = NULL
for (k in 1:200) {
a <- as.data.frame(rdirichlet(100,alpha=floor(runif(n=11,min=0,max=10))))
b <- as.data.frame(rdirichlet(100,alpha=floor(runif(n=11,min=1,max=3))))
c <- as.data.frame(rdirichlet(100,alpha=floor(runif(n=11,min=1,max=11))))
d <- as.data.frame(rdirichlet(100,alpha=floor(runif(n=11,min=1,max=1))))
df <- rbind(df,a,b,c,d)
}
df <-df %>%
rowwise(.) %>%
mutate(samplesize=posamsize(c_across(V1:V11),odds.ratio=2.5,power=0.8,alpha=0.05)$n)