set.seed(156)
n <- 100

### Simuliere Daten y, die nur Fehler sind, und zwei Einflussgrößen x1, x2 mit entsprechend beta1 = beta2 = 0 ###
nrep <- 10000
pvecvoll <- selectedmodelaus2 <- selectedmodelaus4 <- pvecsel <- pvecsel1 <- pvecsel2 <- rep(NA, nrep)

# H_0: beta_1 = 0

for (i in 1:nrep){
  y <- rnorm(n)
  x1 <- rnorm(n)
  x2 <- rnorm(n)  
  
  # bis zu 4 Modelle zur Auswahl
  m1 <- lm(y ~ x1 + x2)
  m2 <- lm(y ~ x1)
  m3 <- lm(y ~ x2)
  m4 <- lm(y ~ 1)
  
  ### 1. Test für H_0 im vollen Modell. Der p-Wert ist gleichverteilt auf [0,1] ###  
  pvecvoll[i] <- summary(m1)$coef[2,4]
  
  ### 2. selektiere die Variable ins Modell, die das kleiner AIC ergibt. p-Werte für diese Variable sind nicht mehr gleichverteilt. ###  
  selectedmodelaus2[i] <- which.min(c(AIC(m2), AIC(m3)))
  pvecsel[i] <- c(summary(m2)$coef[2,4],summary(m3)$coef[2,4])[selectedmodelaus2[i]]
  
  # waehle das beste Modell über AIC
  selectedmodelaus4[i] <- which.min(c(AIC(m1), AIC(m2), AIC(m3), AIC(m4))) 
  # 3. (ii) Test für H_0 im selektierten Modell, nur falls x1 selektiert (bedingt auf Selektion ins Modell)
  pvecsel2[i] <- if(selectedmodelaus4[i]==1) 
    summary(m1)$coef[2,4] else if(selectedmodelaus4[i]==2)  
      summary(m2)$coef[2,4] else NA  
  # 3. (i) Test für H_0 im selektierten Modell, zähle Nichtselektion als Nichtablehnung
  pvecsel1[i] <- if(selectedmodelaus4[i]==1) 
    summary(m1)$coef[2,4] else if(selectedmodelaus4[i]==2)  
      summary(m2)$coef[2,4] else 1
}

table(selectedmodelaus2)
table(selectedmodelaus4)

hist(pvecvoll)
mean(pvecvoll < 0.05, na.rm=T) # Fehler 1. Art im vollen Modell ca. 5 Prozent

hist(pvecsel)
mean(pvecsel < 0.05, na.rm=T) # ca. 10 Prozent signifikant bei alpha = 0.05 

hist(pvecsel1)
mean(pvecsel1 < 0.05, na.rm=T) # (hier) ca. 5 Prozent Ablehnungen, wenn man Nichtselektion ins Modell als Nicht-Ablehnung zählt

hist(pvecsel2)
mean(pvecsel2 < 0.05, na.rm=T) # bedingt auf die Selektion ins Modell erhoehte Fehlerrate 1. Art, ca. 30%


###### selective inference ######

library(selectiveInference)

pvecselFS <- rep(NA,nrep)

for (i in 1:nrep){
  
  y <- rnorm(n)
  x1 <- rnorm(n)
  x2 <- rnorm(n)  
  
  # bis zu 4 Modelle zur Auswahl
  m1 <- lm(y ~ x1 + x2)
  m2 <- lm(y ~ x1)
  m3 <- lm(y ~ x2)
  m4 <- lm(y ~ 1)
  
  
  fsfit=fs(x = matrix(c(x1,x2),ncol=2),y = y, #index = c(1,2), 
           maxsteps = 1)
  pvecselFS[i] <- fsInf(fsfit)$pv
  
}

hist(pvecselFS)
mean(pvecselFS < 0.05)