# Simulate case-control study in R

I'm trying to run a code from chapter 4 of this paper where the response $$y$$ (status of disease) takes value 1 for cases and 0 for controls. Also, $$s_0$$ be the index set of each $$x_{ij}$$ that are associated with the disease and $$n_1=n_2=500$$.

1. For $$j \notin s_0, x_{ij}$$ are simulated under Binomial distribution $$Bin(2,p_j)$$ where $$p_j$$ was generated from $$Beta(\alpha=2,\beta=2)$$ independently for each $$j$$ and in each simulation run.

2. For $$j\in s_0$$, $$x_{ij}$$ were generated in the same way as in control group (what does this mean?)

To simulate from this model in R, I need to find the conditional probabilities and I get

$$P(X(s_0)=x(s_0)|Y=1) = P(X(s_0) = x(s_0)|Y=0)\exp(\alpha^* + x(s_0)^T\beta_0)$$ for normalisation constant $$\alpha^*$$ which i incorrectly calculate as Infinity in my code below.. where i solve $$1/27\exp(\alpha^* + x(s_0)^T\beta_0)=1$$

So the conditional probability when $$Y=0$$ is given by $$1/3^m$$ since when $$s_0$$ contains $$m$$ variables, there are $$3^m$$ possible $$x(s_0)$$ distinct vectors of length $$m$$ with $$0,-1,1$$ entries.

i.e. When $$m$$ = 3, there are $$3^3$$ active $$x(s_0)$$, so we can let the conditional probability = $$1/27$$. The specifications for one model are for $$m=3$$, $$P=500$$ (number of features) and $$\beta_0 = (0.5,0.6,0.7)$$. I then need to sample each of these x(s0) with respect to these probabilities and use it in the package glmnet to model such that such that: MODEL=glmnet(xx, y, family="binomial", alpha=0.99, pmax=40). The vectors for the n2 cases are sampled with replacement, any idea how to do this?

n_cases=500
n_controls=500

#For j not in s0
alpha = 2
beta = 2
pj <- rbeta(500,alpha,beta)
xj <- rbinom(500,2,pj)

#For j in s0
m=3
x = as.matrix(expand.grid(c(-1:1),c(-1:1), c(-1:1)))

beta_0 = c(0.5,0.6,0.7)
Prob = 1/27

##Calculate normalisation Constant - NaN

for (i in 1:3^m){
norm_alpha = log(3^m) - log(sum(x[i,]%*%beta_0))
return(Prob*exp(norm_alpha + x[i,]%*%beta_0))
}

#Calculate Conditional Probability to sample x from:
P_new = Prob*exp(norm_alpha + x%*%beta_0)

#GLMNET Model
X<-sample(xj, replace = TRUE, size = n, prob=P_new) #Wrong size!



So my problems are outlined below:

1. Have i calculated normalisation constant correctly as I get NaNs?
2. How do i sample values of x(s0) to use in glmnet?
3. For values of $$j\notin s_0$$, how/when do i use these?
• What is the distribution of $Y$? – Xi'an Mar 25 at 7:25
• @Xi'an It is Binomially distributed it seems – Btzzzz Mar 25 at 7:35
• It doesn't have a meaningful distribution. You are sampling on (what is assumed to be) the outcome. Typically, all cases are used & a sample w/ the same n is found that somehow 'match' (eg, randomly selected patients that were admitted to the hospital on the same day). – gung - Reinstate Monica Mar 26 at 15:45