# When sigma is known
set.seed(0)
interaction.ind <- t(combn(4,2))
X <- matrix(rnorm(50*4,1,0.1), 50, 4)
epl <- rnorm(50,0,0.01)
y <- 1+X[,1]+X[,2]+X[,1]*X[,2] + epl
ABC(X, y, sigma = 0.01, varind = c(1,2,5), interaction.ind = interaction.ind)
# When sigma is not known
full <- Extract(X, varind = c(1:(dim(X)[2]+dim(interaction.ind)[1])), interaction.ind)
sigma <- selectiveInference::estimateSigma(full, y)$sigmahat # Estimate sigma
Run the code above in your browser using DataLab