2
$\begingroup$

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?
$\endgroup$
  • $\begingroup$ What is the distribution of $Y$? $\endgroup$ – Xi'an Mar 25 at 7:25
  • $\begingroup$ @Xi'an It is Binomially distributed it seems $\endgroup$ – Btzzzz Mar 25 at 7:35
  • $\begingroup$ 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). $\endgroup$ – gung - Reinstate Monica Mar 26 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.