####################### ### Selektionseffekte # von Julien P. Irmer ### Probit-Modell: der Selektionsmechanismus --- set.seed(123456) # Vergleichbarkeit Z <- rnorm(10000) # 10^4 std. normalverteilte Zufallsvariablen simulieren h <- hist(Z, breaks=50, plot=FALSE) cuts <- cut(h\$breaks, c(-4,-1,4)) plot(h, col=cuts, freq = F) lines(x = seq(-4,4,0.01), dnorm(seq(-4,4,0.01)), col = "blue", lwd = 3) abline(v = -1, lwd = 3, lty = 3, col = "gold3") h <- hist(Z, breaks=50, plot=FALSE) cuts <- cut(h\$breaks, c(-4,2,4)) plot(h, col=cuts, freq = F) lines(x = seq(-4,4,0.01), dnorm(seq(-4,4,0.01)), col = "blue", lwd = 3) abline(v = 2, lwd = 3, lty = 3, col = "gold3") set.seed(123456) # Vergleichbarkeit Z <- rnorm(10000) # 10^4 std. normalverteilte Zufallsvariablen simulieren pnorm(q = -1, mean = 0, sd = 1) # Theoretische Wahrscheinlichkeit, dass Z <= -1 mean(Z <= -1) # Empirische (beoabachtete) Wahrscheinlichkeit, dass Z <= -1 pnorm(q = 2, mean = 0, sd = 1) # Theoretische Wahrscheinlichkeit, dass Z <= 2 mean(Z <= 2) # Empirische (beoabachtete) Wahrscheinlichkeit, dass Z <= 2 set.seed(1234567) n <- 1000 Z <- rnorm(n = n, mean = 0, sd = 1) X <- rnorm(n = n, mean = 2, sd = 5) # die Verteilung von X ist beliebig # Selektion beta0 <- -2 beta1 <- .5 s <- beta0 + beta1*X + Z > 0 probit_model <- glm(s ~ 1 + X, family = binomial(link = "probit")) summary(probit_model) ### Einfluss der Selektion --- set.seed(1234567) n <- 1000 X <- rnorm(n = n, mean = 2, sd = 2) Z <- rnorm(n = n, mean = 0, sd = 1) # Populationsregression b0 <- 0.5 b1 <- 1.2 r <- 2 sigma <- 2 eps <- rnorm(n = n, mean = 0, sd = sigma) Y <- b0 + b1*X + r*Z + eps # Populationsregression # Selektion beta0 <- -2 beta1 <- .5 s <- beta0 + beta1*X + Z > 0 # Selektionsmechanismus Y_obs <-rep(NA, length(Y)) Y_obs[s == 1] <- Y[ s== 1] # beobachtbares Y reg_obs <- lm(Y_obs ~ X) reg_obs reg_pop <- lm(Y ~ X) reg_pop # Korrelation mit dem Residuum cor(X[s==1], Z[s==1]) cor(X, Z) cor.test(X[s==1], Z[s==1]) cor.test(X, Z) # Grafik plot(X, Y, pch = 16, col = "skyblue", cex = 1.5) points(X, Y_obs, pch = 16, col = "black") abline(reg_pop, col = "blue", lwd = 5) abline(reg_obs, col = "gold3", lwd = 5) legend(x = "bottomright", legend = c("all", "observed", "regression: all", "regression: observed"), col = c("skyblue", "black", "blue", "gold3"), lwd = c(NA, NA, 5, 5), pch = c(16, 16, NA, NA)) library(sampleSelection) # Paket laden heckman <- heckit(selection = s ~ 1 + X, outcome = Y_obs ~ 1 + X) summary(heckman) coef(summary(heckman)) # Bsp Datensatz laden load(url("https://pandar.netlify.app/post/HeckData.rda")) ## Appendix ### Appendix A: Code zu `R` Grafiken {#AppendixA} set.seed(123456) # Vergleichbarkeit Z <- rnorm(10000) # 10^4 std. normalverteilte Zufallsvariablen simulieren h <- hist(Z, breaks=50, plot=FALSE) cuts <- cut(h\$breaks, c(-4,-1,4)) # Grenzen festlegen für Färbung plot(h, col=cuts, freq = F) lines(x = seq(-4,4,0.01), dnorm(seq(-4,4,0.01)), col = "blue", lwd = 3) abline(v = -1, lwd = 3, lty = 3, col = "gold3") h <- hist(Z, breaks=50, plot=FALSE) cuts <- cut(h\$breaks, c(-4,2,4)) # Grenzen festlegen für Färbung plot(h, col=cuts, freq = F) lines(x = seq(-4,4,0.01), dnorm(seq(-4,4,0.01)), col = "blue", lwd = 3) abline(v = 2, lwd = 3, lty = 3, col = "gold3") plot(X, Y, pch = 16, col = "skyblue", cex = 1.5) points(X, Y_obs, pch = 16, col = "black") abline(reg_pop, col = "blue", lwd = 5) abline(reg_obs, col = "gold3", lwd = 5) legend(x = "bottomright", legend = c("all", "observed", "regression: all", "regression: observed"), col = c("skyblue", "black", "blue", "gold3"), lwd = c(NA, NA, 5, 5), pch = c(16, 16, NA, NA)) ### Appendix B: Heckman Modell simulieren {#AppendixB} simulate_heckman <- function(n, beta0, beta1, b0, b1, r, sigma) { X <- rnorm(n = n, mean = 2, sd = 2) Z <- rnorm(n = n, mean = 0, sd = 1) # Populationsregression eps <- rnorm(n = n, mean = 0, sd = sigma) Y <- b0 + b1*X + r*Z + eps # Populationsregression # Selektion s <- beta0 + beta1*X + Z > 0 # Selektionsmechanismus Y_obs <-rep(NA, length(Y)) Y_obs[s == 1] <- Y[s == 1] # beobachtbares Y df <- data.frame("X" = X, "s" = s, "Y_obs" = Y_obs) return(df) } set.seed(404) # Vergleichbarkeit # Daten simulieren data_heckman <- simulate_heckman(n = 10^5, beta0 = -2, beta1 = 0.5, b0 = 0.5, b1 = 1.2, r = 2, sigma = 2) head(data_heckman) # Daten ansehen # Heckman Modell schätzen für große n heckman_model <- heckit(selection = s ~ 1 + X, outcome = Y_obs ~ 1 + X, data = data_heckman) summary(heckman_model)