seçilim yanlılığının gruplar üzerinden karşılaştırması

• n kişilik bir grupta, ortalama nedensel etkiler \(\mathbb{E}[Y_{1i} - Y_{0i}]\) olarak yazılır, \(\mathbb{E}\) , beklenen (expectation) olarak okunabilir ve ortalama olarak adlandırılabilir.

• Burada ortalama alma her zamanki şekilde yapılır (yani, bireysel sonuçları toplar ve n’e böleriz)

• Sigortanın ortalama nedensel etkisinin araştırılması doğal olarak sigortalı ve sigortasız insan gruplarının ortalama sağlıklarının karşılaştırılmasıyla başlar.

• Bu karşılaştırma, sigorta durumunu belirtmek için 0 ve 1 değerlerini alan bir kukla değişken olan \(D_i\)’nin oluşturulmasıyla kolaylaştırılır.

• Artık sigortalılar arasındaki ortalama için \(\mathbb{E}[Y_{i} | D_i = 1]\) ve sigortasızlar arasındaki ortalama için \(\mathbb{E}[Y_{i} | D_i = 0]\) yazabiliriz. Bu nicelikler sigorta durumuna bağlı ortalamalardır.

• Sigortalılar için ortalama \(Y_i\), zorunlu olarak \(Y_{1i}\) sonucunun bir ortalamasıdır, ancak \(Y_{0i}\) hakkında hiçbir bilgi içermez. Benzer şekilde, sigortasızlar arasındaki ortalama \(Y_i\), \(Y_{0i}\) sonucunun bir ortalamasıdır, ancak bu ortalama, karşılık gelen \(Y_{1i}\) hakkında bilgi içermez.

• Başka bir deyişle, sigortalı olanların izlediği yol \(Y_{1i}\) ile sonlanırken, sigortasız olanların izlediği yol \(Y_{0i}\)’ye çıkar. Bu da sigorta durumuna göre ortalama sağlık farkı hakkında basit ama önemli bir sonuca götürür:

\[ \textrm{Grup Ortalamalarının Farkı} = \mathbb{E}[Y_{i} | D_i = 1] - \mathbb{E}[Y_{i} | D_i = 0] = \mathbb{E}[Y_{1i} | D_i = 1] - \mathbb{E}[Y_{0i} | D_i = 0] \]

• Herkesin \(Y_{1i}\)’sini ve herkesin \(Y_{0i}\)’sini içeren ortalama bir nedensel etki olan \(\mathbb{E}[Y_{1i} - Y_{0i}]\)’nin peşindeyiz, ancak ortalama \(Y_{1i}\)’yi yalnızca sigortalılar için ve ortalama \(Y_{0i}\)’yi yalnızca sigortasızlar için görüyoruz.

• Bu denklemi daha iyi anlamak için, sağlık sigortasının herkesi sabit bir miktarda, \(\kappa\), daha sağlıklı hale getirdiğini hayal etmek yardımcı olur.

\[ Y_{1i} = Y_{0i} + \kappa \]

• Başka bir deyişle, \(\kappa\) sigortanın sağlık üzerindeki hem bireysel hem de ortalama nedensel etkisini ifade eder.

\[ \textrm{Grup Ortalamalarının Farkı} = \mathbb{E}[Y_{0i} + \kappa | D_i = 1] - \mathbb{E}[Y_{0i} | D_i = 0] = \kappa + \mathbb{E}[Y_{0i} | D_i = 1] - \mathbb{E}[Y_{0i} | D_i = 0] \]

• Bu denklem, sigortası olanlar ve olmayanlar arasındaki sağlık farkının; nedensel etki (\(\kappa\)) artı sigortalı ve sigortasız arasındaki ortalama \(Y_{0i}\) farkına eşit olduğunu ortaya koymaktadır.

• Khuzdar ve Maria benzetmesinde olduğu gibi, bu ikinci terim seçim yanlılığını tanımlar. Özellikle, sigorta durumuna göre ortalama sağlık farkı şu şekilde yazılabilir:

\[ \textrm{Grup Ortalamalarının Farkı} = \textrm{Ortalama Nedensel Etki} + \textrm{Seçim Yanlılığı} \]

• Burada seçim yanlılığı, karşılaştırılan gruplar arasındaki ortalama \(Y_{0i}\) farkı olarak tanımlanır.

• Örneğin, sigorta karşılaştırmasında seçim yanlılığının tek kaynağının eğitim olduğunu varsayalım. Bu yanlılık, aynı eğitime sahip kişilerden oluşan örneklere, örneğin üniversite mezunlarına odaklanılarak ortadan kaldırılır.

• Eğitim, böyle bir örneklemdeki sigortalı ve sigortasız kişiler için aynıdır, çünkü örneklemdeki herkes için aynıdır.

• Sağlık sigortası olan ve olmayan kişilerin birçok görünür şekilde farklılık göstermesi, gözlemlenen özellikleri sabit tutsak bile, sigortasızların sigortalılardan muhtemelen bizim görmediğimiz şekillerde farklılık göstereceğini düşündürmektedir (sonuçta, görebildiğimiz değişkenlerin listesi kısmen tesadüfidir).

• Başka bir deyişle, aynı eğitim, gelir ve istihdam durumuna sahip sigortalı ve sigortasız kişilerden oluşan bir örneklemde bile, sigortalıların \(Y_{0i}\) değerleri daha yüksek olabilir.

• Uzmanlarının karşılaştığı temel zorluk, bu tür gözlemlenmemiş farklılıklardan kaynaklanan seçim yanlılığının ortadan kaldırılmasıdır.

Seçim Yanlılığına Sahip Ortalama Farklara Örnek

Bu seçim yanlılığına sahip ortalama farkları bulmak için, ABD nüfusunun yıllık bir anketi olan Ulusal Sağlık Görüşme Anketi’dir (NHIS)’i kullanacağız. Bu anket, Sağlık ve sağlık sigortası hakkında ayrıntılı bilgiler sunar.

• NHIS’de şu soru da vardır: Genel olarak sağlığınızın “mükemmel, çok iyi, iyi, orta veya kötü olduğunu söyler misiniz?

• Bu soru ankette, sigortalı olan veya olmayan, evli 2009 NHIS katılımcılarından oluşan bir örneklemde mükemmel sağlığa 5 ve kötü sağlığa 1 atayan cevaplardan oluşur.

• Bu bağlamda, sigortası olanlar tedavi grubu; sigortası olmayanlar karşılaştırma veya kontrol grubu olarak düşünülebilir.

• Aşağıda ki Tablo 1’in ilk satırı, sigortalı ve sigortasız Amerikalıların ortalama sağlık endeksini karşılaştırır ve istatistikler, eşler için ayrı ayrı tablo halinde sunulur

Tablo 1’in oluşturulması.

Mastering Metrics; Angrist & Pischke Table 1.1

nhis <- read_csv("nhis.csv")

Örnek olarak size bu tabloda bulunan Health Index satırının sonuçlarını çıkarmaya çalışacağım.

• nhis veri setinde bulunan, fml değişkeni, birey kadınsa 1 değerini alan kukla değişkendir.

• hi değişkeni, birey sağlık sigortasına sahipse 1 değerini alan kukla değişkendir.

• hlth değişkeni, bireyin 1’den 5’e kadar (5 en yüksek) sağlık puanını gösterir.

nhis %>% group_by(fml, hi) %>%
  summarize(mean_health = mean(hlth, na.rm=T),
            sd_health = sd(hlth, na.rm=T))

## `summarise()` has grouped output by 'fml'. You can override using the `.groups`
## argument.

## # A tibble: 4 × 4
## # Groups:   fml [2]
##     fml    hi mean_health sd_health
##   <dbl> <dbl>       <dbl>     <dbl>
## 1     0     0        3.70     1.01 
## 2     0     1        3.98     0.934
## 3     1     0        3.61     1.02 
## 4     1     1        3.99     0.928

• Bu ortalamalar tablo ortalamalarıyla örtüşmüyor, çünkü her bir gözlemin kendi özel ağırlığı mevcut, bu ağırlık veri setinde perweight olarak geçiyor.

• weighted.mean fonksiyonu tablodaki sonuçları bize sunuyor.

nhis %>%
  group_by(hi, fml) %>%
  summarise(
    mean_hlth = weighted.mean(hlth, perweight, na.rm = TRUE))

## `summarise()` has grouped output by 'hi'. You can override using the `.groups`
## argument.

## # A tibble: 4 × 3
## # Groups:   hi [2]
##      hi   fml mean_hlth
##   <dbl> <dbl>     <dbl>
## 1     0     0      3.70
## 2     0     1      3.62
## 3     1     0      4.01
## 4     1     1      4.02

Ağırlıklı ortalama nasıl hesaplanır? (weighted.mean fonksiyonu aslında nasıl çalışıyor?)

• Her bir değeri ağırlığıyla çarpın.

nhis$weighted_hlth <- nhis$hlth * nhis$perweight

• Her bir grup için toplam değeri ve toplam ağırlıkları bulun.

sum_weighted_values <- nhis %>% group_by(fml, hi) %>% summarize(toplam_deger = sum(weighted_hlth),
                                                                toplam_agırlık = sum(perweight))

## `summarise()` has grouped output by 'fml'. You can override using the `.groups`
## argument.

sum_weighted_values

## # A tibble: 4 × 4
## # Groups:   fml [2]
##     fml    hi toplam_deger toplam_agırlık
##   <dbl> <dbl>        <dbl>          <dbl>
## 1     0     0     17198932        4653826
## 2     0     1    118121169       29464737
## 3     1     0     14185752        3915407
## 4     1     1    114889967       28598239

• Her bir grup için ortalamayı bulun

agırlıklı_ortalamalar <- sum_weighted_values %>% mutate(ortalamalar =toplam_deger/toplam_agırlık)
agırlıklı_ortalamalar

## # A tibble: 4 × 5
## # Groups:   fml [2]
##     fml    hi toplam_deger toplam_agırlık ortalamalar
##   <dbl> <dbl>        <dbl>          <dbl>       <dbl>
## 1     0     0     17198932        4653826        3.70
## 2     0     1    118121169       29464737        4.01
## 3     1     0     14185752        3915407        3.62
## 4     1     1    114889967       28598239        4.02

Ağırlıklı ortalama bulmak için diğer bir yöntem: Survey paketi

nhis_filtered <- svydesign(ids = ~1, data = nhis, weights = ~perweight)

# Group means and SDs
group_means <- svyby(~hlth, ~hi + fml, nhis_filtered, svymean, na.rm = TRUE) %>%
  rename(mean_hlth = hlth)
group_means

##     hi fml mean_hlth         se
## 0.0  0   0  3.695654 0.03168334
## 1.0  1   0  4.008899 0.01243686
## 0.1  0   1  3.623059 0.03301355
## 1.1  1   1  4.017379 0.01227025

Tablo’nun Health index satırı için, 1., 2., 4. ve 5. sütunlarında bulunan, ortalama değerleri bulduk. Sıra köşeli parantez içindeki standart hataları bulmakta.

Ağırlıklı Standart Hata

nhis_weighted <- nhis %>% group_by(fml, hi) %>% 
  mutate(toplam_deger = sum(weighted_hlth),
         toplam_agırlık = sum(perweight),
         ortalama = toplam_deger/toplam_agırlık) %>% ungroup()

nhis_weighted$sq_dev <- (nhis_weighted$hlth - nhis_weighted$ortalama)^2

nhis_weighted$weighted_sq_dev <- nhis_weighted$sq_dev * (nhis_weighted$perweight / nhis_weighted$toplam_agırlık)

nhis_weighted %>% group_by(fml, hi) %>% summarize(sumweighted_sq_dev = sum(weighted_sq_dev),
                                                  n_effective = sum(perweight)^2 / sum(perweight^2),
                                                  correction_factor = n_effective / (n_effective - 1),
                                                  manual_w_var = sumweighted_sq_dev * correction_factor,
                                                  manual_w_sd = sqrt(manual_w_var))

## `summarise()` has grouped output by 'fml'. You can override using the `.groups`
## argument.

## # A tibble: 4 × 7
## # Groups:   fml [2]
##     fml    hi sumweighted_sq_dev n_effective correction_factor manual_w_var
##   <dbl> <dbl>              <dbl>       <dbl>             <dbl>        <dbl>
## 1     0     0              1.02        1022.              1.00        1.02 
## 2     0     1              0.863       5606.              1.00        0.863
## 3     1     0              1.03         962.              1.00        1.03 
## 4     1     1              0.845       5613.              1.00        0.846
## # ℹ 1 more variable: manual_w_sd <dbl>

Ağırlıklı standart hata: Survey paketi

group_vars <- svyby(~hlth, ~hi + fml, nhis_filtered, svyvar, na.rm = TRUE) %>%
  mutate(sd_health = sqrt(hlth)) %>%
  select(hi, fml, sd_health)
group_vars

##     hi fml sd_health
## 0.0  0   0 1.0121241
## 1.0  1   0 0.9288020
## 0.1  0   1 1.0140217
## 1.1  1   1 0.9195026

group_summary <- left_join(group_means, group_vars, by = c("hi", "fml"))
group_summary

##   hi fml mean_hlth         se sd_health
## 1  0   0  3.695654 0.03168334 1.0121241
## 2  1   0  4.008899 0.01243686 0.9288020
## 3  0   1  3.623059 0.03301355 1.0140217
## 4  1   1  4.017379 0.01227025 0.9195026

3. ve 6. sütunlardaki, Health index farkları (difference)

Farklar sadece ortalama farklar bulunarak gösterilebilir

wife_diff <- 4.017379    - 3.623059 
husbands_diff <- 4.008899 - 3.695654
wife_diff; husbands_diff

## [1] 0.39432

## [1] 0.313245

Farkların diğer gösterimi Regresyon yöntemi

wife_lm <- lm(hlth ~ hi, data = nhis %>% filter(fml == TRUE), weights = perweight)
husband_lm <- lm(hlth ~ hi, data = nhis %>% filter(fml == FALSE), weights = perweight)
summary(wife_lm); summary(husband_lm)

## 
## Call:
## lm(formula = hlth ~ hi, data = nhis %>% filter(fml == TRUE), 
##     weights = perweight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -409.48  -38.47   -0.91   47.57  178.60 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.62306    0.02769  130.84   <2e-16 ***
## hi           0.39432    0.02953   13.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.79 on 9393 degrees of freedom
## Multiple R-squared:  0.01864,    Adjusted R-squared:  0.01853 
## F-statistic: 178.4 on 1 and 9393 DF,  p-value: < 2.2e-16

## 
## Call:
## lm(formula = hlth ~ hi, data = nhis %>% filter(fml == FALSE), 
##     weights = perweight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -347.42  -41.47   -0.47   49.69  175.89 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.69565    0.02628  140.65   <2e-16 ***
## hi           0.31325    0.02827   11.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 56.68 on 9393 degrees of freedom
## Multiple R-squared:  0.0129, Adjusted R-squared:  0.01279 
## F-statistic: 122.7 on 1 and 9393 DF,  p-value: < 2.2e-16

• hi katsayısı farkı (Difference) sütunlarını verir. (Intercept) katsayısı, yukarıda bulduğumuz, sigorta olmayan NO HI sütunlarını verir.

• hi katsayısı ve (Intercept) katsayısının toplamı, **Some HI* sütunlarını verir.

Survey paketi kullanarak Diffrence sütununu ve standart hatasını bulmak

# svyttest for differences
ttest_women <- svyttest(hlth ~ hi, subset(nhis_filtered, fml == TRUE))
ttest_men <- svyttest(hlth ~ hi, subset(nhis_filtered, fml == FALSE))
ttest_women; ttest_men

## 
##  Design-based t-test
## 
## data:  hlth ~ hi
## t = 11.196, df = 9393, p-value < 2.2e-16
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
##  0.3252807 0.4633587
## sample estimates:
## difference in mean 
##          0.3943197

## 
##  Design-based t-test
## 
## data:  hlth ~ hi
## t = 9.2031, df = 9393, p-value < 2.2e-16
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
##  0.2465255 0.3799649
## sample estimates:
## difference in mean 
##          0.3132452

• Difference in mean satırının altında, grup için sigortalı olmanın ve olmamanın farkı bulunur (difference sütunları).

• Difference in mean tahmini t değerine bölersek, standart hataya ulaşırız.

# Extracting scalars (ensure they are not NULL)
sd_diff_wife <- as.numeric(ttest_women$estimate/ttest_women$statistic)
sd_diff_husband <- as.numeric(ttest_men$estimate/ttest_men$statistic)
sd_diff_wife; sd_diff_husband

## [1] 0.03522008

## [1] 0.03403689

diff_wife <- as.numeric(ttest_women$estimate)
diff_husband <- as.numeric(ttest_men$estimate)
diff_wife; diff_husband

## [1] 0.3943197

## [1] 0.3132452

Benzer bir regresyonu Survey paketiyle de yapabiliriz.

summary(svyglm(hlth ~ hi, subset(nhis_filtered, fml == TRUE)))

## 
## Call:
## svyglm(formula = hlth ~ hi, design = subset(nhis_filtered, fml == 
##     TRUE))
## 
## Survey design:
## subset(nhis_filtered, fml == TRUE)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.62306    0.03301   109.7   <2e-16 ***
## hi           0.39432    0.03522    11.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.8674062)
## 
## Number of Fisher Scoring iterations: 2

summary(svyglm(hlth ~ hi, subset(nhis_filtered, fml == FALSE)))

## 
## Call:
## svyglm(formula = hlth ~ hi, design = subset(nhis_filtered, fml == 
##     FALSE))
## 
## Survey design:
## subset(nhis_filtered, fml == FALSE)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.69565    0.03168 116.643   <2e-16 ***
## hi           0.31325    0.03404   9.203   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.8846404)
## 
## Number of Fisher Scoring iterations: 2

• Bu durumda svyttest’e ihtiyaç kalmayabilir, çünkü doğru standart hatalara bu yöntemle de ulaşılabilir.

# First, run the regressions for males and females
male_reg <- svyglm(hlth ~ hi, design = subset(nhis_filtered, fml == FALSE))
female_reg <- svyglm(hlth ~ hi, design = subset(nhis_filtered, fml == TRUE))

# Extract coefficients and standard errors
male_coef <- tidy(male_reg)
female_coef <- tidy(female_reg)

# Function to calculate means and SDs for specific groups
get_sd <- function(design, var_name, hi_val, fml_val) {
  group_design <- subset(design, hi == hi_val & fml == fml_val)
  sd_val <- sqrt(svyvar(as.formula(paste0("~", var_name)), group_design))[1]
  return(sd_val)
}

# Get standard deviations for each group
male_hi0_sd <- get_sd(nhis_filtered, "hlth", 0, FALSE)
male_hi1_sd <- get_sd(nhis_filtered, "hlth", 1, FALSE)
female_hi0_sd <- get_sd(nhis_filtered, "hlth", 0, TRUE)
female_hi1_sd <- get_sd(nhis_filtered, "hlth", 1, TRUE)

# Create a data frame for the first row (health variable)
health_row <- data.frame(
  Variable = "hlth",
  
  # Males (fml = 0)
  Male_hi1 = paste0(round(male_coef$estimate[1] + male_coef$estimate[2], 2), 
                   " [", round(male_hi1_sd, 2), "]"),
  Male_hi0 = paste0(round(male_coef$estimate[1], 2), 
                   " [", round(male_hi0_sd, 2), "]"),
  Male_diff = paste0(round(male_coef$estimate[2], 2), 
                    " (", round(male_coef$std.error[2], 2), ")"),
  
  # Females (fml = 1)
  Female_hi1 = paste0(round(female_coef$estimate[1] + female_coef$estimate[2], 2), 
                     " [", round(female_hi1_sd, 2), "]"),
  Female_hi0 = paste0(round(female_coef$estimate[1], 2), 
                     " [", round(female_hi0_sd, 2), "]"),
  Female_diff = paste0(round(female_coef$estimate[2], 2), 
                      " (", round(female_coef$std.error[2], 2), ")")
)

# Create the gt table
table <- gt(health_row) %>%
  tab_header(
    title = "Health Outcomes by Insurance Status and Gender"
  ) %>%
  cols_label(
    Variable = "",
    Male_hi1 = "Some HI",
    Male_hi0 = "No HI",
    Male_diff = "Difference",
    Female_hi1 = "Some HI",
    Female_hi0 = "No HI",
    Female_diff = "Difference"
  ) %>%
  tab_spanner(
    label = "Husbands",
    columns = c(Male_hi1, Male_hi0, Male_diff)
  ) %>%
  tab_spanner(
    label = "Wives",
    columns = c(Female_hi1, Female_hi0, Female_diff)
  ) %>%
  tab_footnote(
    footnote = "Standard deviations in square brackets; Standard errors in parentheses.",
    locations = cells_title()
  )

# Display the table
table

# Function to add more variables to the table
add_variable_row <- function(table_data, var_name, nhis_filtered, digit=2) {
  # Run regressions
  male_reg <- svyglm(as.formula(paste0(var_name, " ~ hi")), 
                     design = subset(nhis_filtered, fml == FALSE))
  female_reg <- svyglm(as.formula(paste0(var_name, " ~ hi")), 
                       design = subset(nhis_filtered, fml == TRUE))
  
  # Extract coefficients
  male_coef <- tidy(male_reg)
  female_coef <- tidy(female_reg)
  
  # Get standard deviations
  male_hi0_sd <- get_sd(nhis_filtered, var_name, 0, FALSE)
  male_hi1_sd <- get_sd(nhis_filtered, var_name, 1, FALSE)
  female_hi0_sd <- get_sd(nhis_filtered, var_name, 0, TRUE)
  female_hi1_sd <- get_sd(nhis_filtered, var_name, 1, TRUE)
  
  # Create a new row
  new_row <- data.frame(
    Variable = var_name,
    Male_hi1 = paste0(round(male_coef$estimate[1] + male_coef$estimate[2], digit), 
                     " [", round(male_hi1_sd, digit), "]"),
    Male_hi0 = paste0(round(male_coef$estimate[1], digit), 
                     " [", round(male_hi0_sd, digit), "]"),
    Male_diff = paste0(round(male_coef$estimate[2], digit), 
                      " (", round(male_coef$std.error[2], digit), ")"),
    Female_hi1 = paste0(round(female_coef$estimate[1] + female_coef$estimate[2], digit), 
                       " [", round(female_hi1_sd, digit), "]"),
    Female_hi0 = paste0(round(female_coef$estimate[1], digit), 
                       " [", round(female_hi0_sd, digit), "]"),
    Female_diff = paste0(round(female_coef$estimate[2], digit), 
                        " (", round(female_coef$std.error[2], digit), ")")
  )
  
  # Bind to existing data
  rbind(table_data, new_row)
}

# To add x and z variables:
# all_rows <- health_row
# all_rows <- add_variable_row(all_rows, "x", nhis_filtered)
# all_rows <- add_variable_row(all_rows, "z", nhis_filtered)
# 
# full_table <- gt(all_rows) %>%
#   # (Same formatting as above)

varlist <- c(“hlth”, “nwhite”, “age”, “yedu”, “famsize”, “empl”, “inc”)

all_rows <- health_row
all_rows <- add_variable_row(all_rows, "nwhite", nhis_filtered)
all_rows <- add_variable_row(all_rows, "age", nhis_filtered)
all_rows <- add_variable_row(all_rows, "yedu", nhis_filtered)
all_rows <- add_variable_row(all_rows, "famsize", nhis_filtered)
all_rows <- add_variable_row(all_rows, "empl", nhis_filtered)
all_rows <- add_variable_row(all_rows, "inc", nhis_filtered, digit=0)

all_rows$Variable <- c("Health Index", "Nonwhite", "Age", "Education", "Family size", "Employed", "Family income")

full_table <- gt(all_rows)  %>%
  tab_header(
    title = "Health Outcomes by Insurance Status and Gender"
  ) %>%
  cols_label(
    Variable = "",
    Male_hi1 = "Some HI",
    Male_hi0 = "No HI",
    Male_diff = "Difference",
    Female_hi1 = "Some HI",
    Female_hi0 = "No HI",
    Female_diff = "Difference"
  ) %>%
  tab_spanner(
    label = "Husbands",
    columns = c(Male_hi1, Male_hi0, Male_diff)
  ) %>%
  tab_spanner(
    label = "Wives",
    columns = c(Female_hi1, Female_hi0, Female_diff)
  ) %>%
  tab_footnote(
    footnote = "Standard deviations in square brackets; Standard errors in parentheses.",
    locations = cells_title()
  )

full_table

• Tablodan şunu görüyoruz: Sigortalılar daha yaşlı, daha eğitimli ve daha zengin. Sağlık farkı sigortadan mı, yoksa bu özelliklerden mi geliyor?

Peki bu seçim yanlılığını nasıl yeneriz? Tabii ki rastgele atama ile! Bir deneyde, kimin sigorta alacağına yazı tura ile karar verirsek, gruplar birbirine benzer hale gelir. Farklılıklar sadece sigortadan kaynaklanır. Mesela, bir partide kimin kahve içeceğini rastgele seçsek, konuşkanlığın kahveden mi yoksa karakterden mi geldiğini anlarız.

RAND Sağlık Sigortası Deneyi

Peki, bu yazı tura fikri gerçek hayatta işe yaradı mı? Evet, hem de nasıl! 1974-1982 yılları arasında Amerika’da RAND Sağlık Sigortası Deneyi (HIE) yapıldı—ekonometrinin Yıldız Savaşları destanı gibi bir şey! 3.958 kişi, ülkenin farklı yerlerinden toplandı ve 14 farklı sigorta planına rastgele atandı. Bazıları her şeyi bedava aldı (ücretsiz plan, Free Plan), bazıları ise “Paran yoksa doktora gitme!” tarzı felaket planlarına (Catastrophic Plan) düştü. Arada da muafiyetli (Deductible) ve eş sigortalı (Coinsurance) planlar vardı—herkesin bütçesine göre bir şeyler!

HIE’nin amacı neydi? Ekonomistler şunu merak ediyordu: “Sağlık hizmeti ucuzladıkça insanlar doktora koşar mı?” ve “Daha iyi sigorta gerçekten daha sağlıklı yapar mı?” İlk sorunun cevabı kocaman bir “EVET!”—ücretsiz plan alanlar doktora adeta maraton koşar gibi gitti. Ama ikinci soru? Orası biraz bulanık. Sağlık sonuçları o kadar da değişmedi. Peki, bu nasıl oldu? Hadi verilere bakalım!

HIE’ye ülkenin altı bölgesinden 14 ila 61 yaş arası 3.958 kişi kaydoldu. HIE örneği Medicare katılımcılarını ve çoğu Medicaid ve askeri sağlık sigortası abonelerini çalışmadan hariç tuttu. HIE katılımcıları 14 sigorta planından birine rastgele atandı.

Kitapda (Mastering Metrics) 14 orjinal plan yerine dört plan grubuyla çalışılmış: felaket Catastrophic Plan, muafiyetli Deductible Plan, eş sigortalı Coinsurance Plan ve ücretsiz Free Plan.

Felaket planları, maliyetin neredeyse tamamını sigorta sahibine yüklerken, muafiyet, eş ödeme ve ücretsiz planlar artan kapsam seviyeleriyle karakterize edilir.

HIE Örneği

Şimdi HIE verilerine geçelim—burada sigorta rastgele dağıtıldı, yani seçim yanlılığı yok!

Rastgele dağıtıldı dedik ama, bundan emin miyiz? Bu dört sağlık sigortası planının, karşılaştırılabilir kişilere dağıtılıp dağıtılmadığını anlamak için, denge kontrolü yaparız—grupların özelliklerini karşılaştırırız.

Mastering Metrics; Angrist & Pischke Table 1.3

Bu tabloyu oluşturmak istiyoruz.

rand_clean <- read_csv("rand_clean.csv")

table(rand_clean$pt)

## 
## Catastrophic  Coinsurance   Deductible         Free 
##          759         1022          881         1295

rand_clean %>% group_by(pt) %>% summarize(mean_female = mean(female, na.rm = T),
                                          sd_female = sd(female, na.rm = T))

## # A tibble: 4 × 3
##   pt           mean_female sd_female
##   <chr>              <dbl>     <dbl>
## 1 Catastrophic       0.560     0.497
## 2 Coinsurance        0.535     0.499
## 3 Deductible         0.537     0.499
## 4 Free               0.522     0.500

Tablonun ilk satırı, Female, Catastrophic Plan için ortalamayı gösteriyor. Diğer sigorta planları için ise bu ortalamadan ne kadar farklı olduklarını. Catastrophic Plan’da bulunan kadın ortalaması %56’mış, Coinsurance Plan ise Catastrophic Plan’dan 0.025 daha düşük bir kadın ortalamasına sahip.

0.5599473 - 0.5352250

## [1] 0.0247223

Diğer planlar arasındaki farklar da aynı şekilde bulunabilir.

Peki parantez içinde bulunan standart hataları nasıl bulacağız ve ne anlama geliyorlar.

lm_result <- lm(female ~ pt, data = rand_clean)
summary(lm_result)

## 
## Call:
## lm(formula = female ~ pt, data = rand_clean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5599 -0.5352  0.4400  0.4648  0.4780 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.55995    0.01810  30.929   <2e-16 ***
## ptCoinsurance -0.02472    0.02390  -1.034   0.3010    
## ptDeductible  -0.02306    0.02470  -0.933   0.3506    
## ptFree        -0.03794    0.02280  -1.664   0.0962 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4988 on 3953 degrees of freedom
## Multiple R-squared:  0.0007013,  Adjusted R-squared:  -5.713e-05 
## F-statistic: 0.9247 on 3 and 3953 DF,  p-value: 0.4279

Eğer regresyon yaparsak, (Intercept) bize Catastrophic Plan’nın ortalamasını veriyor, Std. Error ise parantez içinde bulunan standart hataları. Catastrophic Plan’nın diğer planlardan farkı ve standart hataları, diğer regresyon katsayılarında bulunabilir.

Plan tiplerinin bulunduğu, pt değişkeninin, referans gözlemini “Catastrophic” yaptığımı gözden kaçırmayın. **rand_clean$pt <- relevel(rand_clean$pt, ref = "Catastrophic")**

Bulmamız gereken bir sütun daha kaldı. O da Catastrophic Plan ile diğer bütün sigorta planlarının birlikte olduğu durumun karşılaştırılması.

Bunu gerçekleştirmek için pt2 kukla değişkeni oluşturuyoruz. Bu kukla değişkenin gözlemleri, pt değişkeni Catastrophic’e eşitse, Catastrophic olacak, değilse Other olarak adlandırılacak.

rand_clean$pt2 <- ifelse(rand_clean$pt == "Catastrophic", "Catastrophic", "Other")

# Convert to factor
rand_clean$pt2 <- factor(rand_clean$pt2)

lm(female ~ pt2, data = rand_clean) %>% summary()

## 
## Call:
## lm(formula = female ~ pt2, data = rand_clean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5599 -0.5303  0.4400  0.4697  0.4697 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.55995    0.01810  30.934   <2e-16 ***
## pt2Other    -0.02962    0.02014  -1.471    0.141    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4987 on 3955 degrees of freedom
## Multiple R-squared:  0.0005467,  Adjusted R-squared:  0.000294 
## F-statistic: 2.163 on 1 and 3955 DF,  p-value: 0.1414

Gördüğünüz gibi son sütunun sonuçlarını da bulduk.

# Define variable list
varlist <- c("female", "blackhisp", "age", "educper",
             "income1cpi", "hosp", "ghindx", "cholest", "diastol",
             "systol", "mhi")

# Run regressions for pt (Catastrophic vs others) and pt2 (All vs Catastrophic)
results <- map_dfr(varlist, function(var) {
  # For pt (Catastrophic vs other categories)
  formula_pt <- as.formula(paste(var, "~ pt"))
  model_pt <- lm(formula_pt, data = rand_clean)
  tidy_pt <- tidy(model_pt) %>%
    mutate(term = recode(term,
                         `(Intercept)` = "Catastrophic",
                         `ptFree` = "Free vs Catastrophic",
                         `ptDeductible` = "Deductible vs Catastrophic",
                         `ptCoinsurance` = "Coinsurance vs Catastrophic"),
           outcome = var)
  
  # For pt2 (All vs Catastrophic)
  formula_pt2 <- as.formula(paste(var, "~ pt2"))
  model_pt2 <- lm(formula_pt2, data = rand_clean)
  tidy_pt2 <- tidy(model_pt2) %>%
    mutate(term = recode(term,
                         `(Intercept)` = "Catastrophic",
                         `pt2Other` = "All vs Catastrophic"),
           outcome = var) %>%
    filter(term != "Catastrophic")  # Remove duplicate "Catastrophic" row from pt2
  
  # Combine both results
  bind_rows(tidy_pt, tidy_pt2)
})

# Relabel terms and add significance stars for both models
results_clean <- results %>%
  mutate(stars = case_when(
    p.value < 0.001 ~ "***",
    p.value < 0.01  ~ "**",
    p.value < 0.05  ~ "*",
    TRUE            ~ ""),
         est_se = sprintf("%.3f%s\n(%.3f)", estimate, stars, std.error)) %>%
  select(outcome, term, est_se)

# Reshape to wide format
results_wide <- results_clean %>%
  pivot_wider(names_from = term, values_from = est_se) %>%
  arrange(outcome)

# Format with gt
results_wide %>%
  gt(rowname_col = "outcome") %>%
  tab_header(
    title = "Regression Results: Variables Regressed on Plan Type",
    subtitle = "Estimates with standard errors in parentheses. * p<0.05, ** p<0.01, *** p<0.001"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(data_row.padding = px(3))

Ortalamalardaki bir farkın standart hatası, istatistiksel kesinliğinin bir ölçüsüdür: Örnek ortalamalarındaki bir fark yaklaşık iki standart hatadan daha küçük olduğunda, farkın genellikle bu örneklerin çekildiği popülasyonların aslında aynı olduğu hipoteziyle uyumlu bir şans bulgusu olduğu düşünülür.

Yaklaşık iki standart hatadan daha büyük olan farkların istatistiksel olarak anlamlı olduğu söylenir: bu gibi durumlarda, bu farkların tamamen şans eseri ortaya çıkması çok olası değildir (ancak imkansız değildir).

İstatistiksel anlamlılık kavramı, Yukarıda bulunan Tablo’daki gibi karşılaştırmaları yorumlamamıza yardımcı olur.

Bu tablodaki farklar çoğunlukla küçüktür.

Tablo’nun son 5 satırı, tedavi grupları arasında tedavi öncesi sonuçlarda makul derecede iyi bir denge olduğuna dair kanıtlarla da sunar.

Bu panel, genel sağlık için tedavi öncesi endekste istatistiksel olarak anlamlı bir fark göstermez.

Benzer şekilde, tedavi öncesi kolesterol, kan basıncı ve ruh sağlığı, istatistiksel anlamlılığa yakın sadece birkaç karşıtlıkla tedavi atamasıyla büyük ölçüde ilgisiz görünmektedir.

Ek olarak, serbest gruptaki düşük kolesterol, felaket grubundakinden biraz daha iyi bir sağlık olduğunu gösterse de, bu iki grup arasındaki genel sağlık endeksindeki farklılıklar ters yöndedir (çünkü daha düşük endeks değerleri daha kötü sağlık olduğunu gösterir).

Tutarlı bir örüntünün olmaması, bu boşlukların şansa bağlı olduğu fikrini güçlendirir.

HIE deneyinin sonuç tahminleri

HIE’den çıkan ilk önemli bulgu, daha cömert sigorta planlarına atanan deneklerin önemli ölçüde daha fazla sağlık hizmeti kullanmasıydı.

Ekonomistlerin, bir malın talebinin ucuzladığında artması gerektiği görüşünü doğrulayan bu bulgu, aşağıdaki tablo’nun A panelinde görülebilir.

Beklenebileceği gibi, muhtemelen kabul kararları genellikle doktorlar tarafından verildiği için, hastanede yatan hasta kabulü ayakta tedaviden daha az fiyata duyarlıdır.

Öte yandan, ücretsiz bakım planına atama, ayakta tedavi harcamalarını felaket planlarındaki harcamalara kıyasla üçte iki oranında (169/248) artırırken, toplam tıbbi harcamalar %45 artar.

Bu büyük farklar ekonomik olarak önemli olduğu kadar istatistiksel olarak da anlamlıdır. Sağlık hizmetlerinin maliyeti konusunda endişelenmek zorunda olmayan denekler açıkça çok daha fazlasını tüketmişlerdir.

Ekstra bakım ve masraf onları daha sağlıklı mı yapıyor?

HIE tedavi grupları arasında sağlık göstergelerini karşılaştıran aşağıdaki Tablo’daki Panel B, cevabın hayır olduğunu öne sürüyor.

Kolesterol seviyeleri, kan basıncı ve genel sağlık ve ruh sağlığının özet endeksleri gruplar arasında dikkat çekici derecede benzerdir (bu sonuçlar çoğunlukla rastgele atamadan 3 veya 5 yıl sonra ölçülmüştür).

Resmi istatistiksel testler, tablodan görülebileceği gibi istatistiksel olarak anlamlı bir fark göstermemektedir.

Bu HIE bulguları birçok ekonomisti cömert sağlık sigortasının istenmeyen sonuçlara yol açabileceğine, daha iyi sağlık şeklinde bir temettü üretmeden sağlık hizmeti kullanımını ve maliyetlerini artırabileceğine ikna etmiştir.

Gelin bu sonuçları biz de bulalım. Deney sonrası sağlık ölçütleri, aynı veri setinde sonlarına x eklenerek kaydedilmişler. Bu sonuçlar, Mastering Metrics Table 1.4’de Panel B olarak sunulmuş.

varlist <- c("ghindxx", "cholestx", "diastolx", "systolx", "mhix")

# Run regressions for pt (Catastrophic vs others) and pt2 (All vs Catastrophic)
results <- map_dfr(varlist, function(var) {
  # For pt (Catastrophic vs other categories)
  formula_pt <- as.formula(paste(var, "~ pt"))
  model_pt <- lm(formula_pt, data = rand_clean)
  tidy_pt <- tidy(model_pt) %>%
    mutate(term = recode(term,
                         `(Intercept)` = "Catastrophic",
                         `ptFree` = "Free vs Catastrophic",
                         `ptDeductible` = "Deductible vs Catastrophic",
                         `ptCoinsurance` = "Coinsurance vs Catastrophic"),
           outcome = var)
  
  # For pt2 (All vs Catastrophic)
  formula_pt2 <- as.formula(paste(var, "~ pt2"))
  model_pt2 <- lm(formula_pt2, data = rand_clean)
  tidy_pt2 <- tidy(model_pt2) %>%
    mutate(term = recode(term,
                         `(Intercept)` = "Catastrophic",
                         `pt2Other` = "All vs Catastrophic"),
           outcome = var) %>%
    filter(term != "Catastrophic")  # Remove duplicate "Catastrophic" row from pt2
  
  # Combine both results
  bind_rows(tidy_pt, tidy_pt2)
})

# Relabel terms and add significance stars for both models
results_clean <- results %>%
  mutate(stars = case_when(
    p.value < 0.001 ~ "***",
    p.value < 0.01  ~ "**",
    p.value < 0.05  ~ "*",
    TRUE            ~ ""),
         est_se = sprintf("%.3f%s\n(%.3f)", estimate, stars, std.error)) %>%
  select(outcome, term, est_se)

# Reshape to wide format
results_wide <- results_clean %>%
  pivot_wider(names_from = term, values_from = est_se) %>%
  arrange(outcome)

# Format with gt
results_wide %>%
  gt(rowname_col = "outcome") %>%
  tab_header(
    title = "Regression Results: Variables Regressed on Plan Type, Panel B",
    subtitle = "Estimates with standard errors in parentheses. * p<0.05, ** p<0.01, *** p<0.001"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(data_row.padding = px(3))

Ancak Panel A için başka veriseti yüklemeliyiz.

combined_df <- read_csv("combined_df.csv")

count(combined_df, plantype)

## # A tibble: 4 × 2
##   plantype         n
##   <chr>        <int>
## 1 Catastrophic  3724
## 2 CostSharing   5464
## 3 Deductible    4175
## 4 Free          6840

varlist <- c("ftf", "out_inf", "totadm", "inpdol_inf", "tot_inf")

# Run regressions for pt (Catastrophic vs others) and pt2 (All vs Catastrophic)
results <- map_dfr(varlist, function(var) {
  # For pt (Catastrophic vs other categories)
  formula_pt <- as.formula(paste(var, "~ plantype"))
  model_pt <- lm(formula_pt, data = combined_df)
  tidy_pt <- tidy(model_pt) %>%
    mutate(term = recode(term,
                         `(Intercept)` = "Catastrophic",
                         `plantypeFree` = "Free vs Catastrophic",
                         `plantypeDeductible` = "Deductible vs Catastrophic",
                         `plantypeCostSharing` = "Coinsurance vs Catastrophic"),
           outcome = var)
  
  # For pt2 (All vs Catastrophic)
  formula_pt2 <- as.formula(paste(var, "~ any_ins"))
  model_pt2 <- lm(formula_pt2, data = combined_df)
  tidy_pt2 <- tidy(model_pt2) %>%
    mutate(term = recode(term,
                         `(Intercept)` = "Catastrophic",
                         `any_insTRUE` = "All vs Catastrophic"),
           outcome = var) %>%
    filter(term != "Catastrophic")  # Remove duplicate "Catastrophic" row from pt2
  
  # Combine both results
  bind_rows(tidy_pt, tidy_pt2)
})

# Relabel terms and add significance stars for both models
results_clean <- results %>%
  mutate(stars = case_when(
    p.value < 0.001 ~ "***",
    p.value < 0.01  ~ "**",
    p.value < 0.05  ~ "*",
    TRUE            ~ ""),
         est_se = sprintf("%.3f%s\n(%.3f)", estimate, stars, std.error)) %>%
  select(outcome, term, est_se)

# Reshape to wide format
results_wide <- results_clean %>%
  pivot_wider(names_from = term, values_from = est_se) %>%
  arrange(outcome)

# Format with gt
results_wide %>%
  gt(rowname_col = "outcome") %>%
  tab_header(
    title = "Regression Results: Variables Regressed on Plan Type, Panel A",
    subtitle = "Estimates with standard errors in parentheses. * p<0.05, ** p<0.01, *** p<0.001"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(data_row.padding = px(3))

Bu tablo bize ne diyor? Farklar genelde küçük ve istatistiksel olarak anlamlı değil—yani rastgelelik işe yaramış! Gruplar, birbiriyle aynı meyve sepetinden çıkmış gibi görünüyor.

HIE’nin sonuçlarına hızlıca bakalım: Ücretsiz plan alanlar sağlık hizmetini çıldırmış gibi kullandı—doktor ziyaretleri, hastane masrafları tavan yaptı! Ama sağlık sonuçları? Kolesterol, kan basıncı, genel sağlık… Hepsi差不多 (neredeyse aynı)! Bu ne anlama geliyor? Daha fazla sigorta, daha fazla doktor ziyareti demek, ama süper sağlıklı bir hayatın garantisi değil.

• Kaynak: Bu notlar, Joshua D. Angrist ve Jörn-Steffen Pischke’nin Mastering Metrics kitabının 1. bölümünden esinlenerek, öğrenciler için eğlenceli bir şekilde uyarlanmıştır.