키보드를 샀다. 노트북 거치대를 사용하면서 노트북 자판을 기울인 채로 타자를 치다보니 손목이 아파서… 큰 맘 먹고 돈을 좀 주고 키보드를 바꿨다. 막상 써보니깐 왜 그동안 나는 바보같이 이 구매를 미뤄왔나 싶다. 어떻게든 이 키보드를 더 두들기고 싶어서 아래 작업을 혼자 뚱땅 뚱땅 해보았다. 도각 도각 소리가 기분이 좋다.

getwd()

## [1] "C:/RRR"

load("C:/RRR/kgss/kgsslogistic_210304.RData")
library(mice)

## 
## Attaching package: 'mice'

## The following object is masked from 'package:stats':
## 
##     filter

## The following objects are masked from 'package:base':
## 
##     cbind, rbind

library(moonBook)

## Warning: package 'moonBook' was built under R version 4.0.4

## 
## Attaching package: 'moonBook'

## The following object is masked from 'package:mice':
## 
##     densityplot

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(ggplot2)
library(naniar)
library(ggeffects)
library(effects)

## Loading required package: carData

## lattice theme set by effectsTheme()
## See ?effectsTheme for details.

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

library(survival)

## Warning: package 'survival' was built under R version 4.0.4

library(broom)
library(descr)

데이터

데이터는 KGSS 08년, 18년 데이터다. 비규범적 성관계에 대한 질문은 08, 13, 18년에 조사되었다. 그런데 13년 데이터를 뺀 이유는 첫째로 13년 데이터는 종교 모듈이 없다. 이 문제는 이번 작업에서는 그다지 중요하진 않다. 사실 더 큰 두 번째 문제는 13년 데이터는 성 역할 규범을 측정할 수 있는 문항이 없다. 그런데 주지하다시피 성관계에 대한 편견, 태도는 성 역할 규범에 아주 큰 영향을 받는다. 그래서 08, 18년 데이터만 사용하였다.

mice 패키지를 사용해서 결측값들을 imputation했다.
관련 설명은 다음 글 참고

Step 1. 변수

종속변수는 각각 혼전 성관계, 혼외 성관계, 동성간 성관계에 대한 4점짜리 척도다. 해당 행위가 얼마나 “잘못되었는가”에 대한 답변이다.

* 혼전 성관계 각각 4첨 척도 빈도

descr::freq(KG4$SEXATT1_n)

## KG4$SEXATT1_n 
##                   Frequency Percent
## 1_Not Wrong             646   25.44
## 2_Sometimes Wrong       827   32.57
## 3_Mostly Wrong          450   17.72
## 4_Totally Wrong         616   24.26
## Total                  2539  100.00

* 혼외 성관계 각각 4첨 척도 빈도

descr::freq(KG4$SEXATT2_n)

## KG4$SEXATT2_n 
##                   Frequency Percent
## 1_Not Wrong              75   2.954
## 2_Sometimes Wrong       176   6.932
## 3_Mostly Wrong          450  17.724
## 4_Totally Wrong        1838  72.391
## Total                  2539 100.000

* 동성간 성관계 각각 4첨 척도 빈도

descr::freq(KG4$SEXATT3_n)

## KG4$SEXATT3_n 
##                   Frequency Percent
## 1_Not Wrong             325   12.80
## 2_Sometimes Wrong       260   10.24
## 3_Mostly Wrong          355   13.98
## 4_Totally Wrong        1599   62.98
## Total                  2539  100.00

설명변수는 1)성별(SEX_n) 2)연령(AGE_n) 3)교육수준(EDUC_n) 4)정치성향(PARTYLR_n) 5)혼인(MARITAL_n) 6)종교(RELIG_n) 6)년도(YEAR_n2) 7)국힘당 지지(party_n1) 8)더민주 지지(party_n2) 9)성역할규범(FAMILY_n2)이다.

# Step 2. 순서형 로짓 분석 ## 혼전 성관계

r reg1 <- polr(SEXATT1_n ~ SEX_n + AGE_n + MARITAL_n + party_n1 + party_n2 + PARTYLR_n + EDUC_n + FAMILY_n2 + YEAR_n2*RELIG_n ,data = KG4, Hess = T)

혼외 성관계

reg2 <- polr(SEXATT2_n ~ SEX_n + AGE_n + MARITAL_n  + party_n1 + party_n2
             + PARTYLR_n + EDUC_n
             + FAMILY_n2 + YEAR_n2*RELIG_n
             ,data = KG4, Hess = T)

동성간 성관계

reg3 <- polr(SEXATT3_n ~ SEX_n + AGE_n + MARITAL_n + party_n1 + party_n2
             + PARTYLR_n + EDUC_n
             + FAMILY_n2 + YEAR_n2*RELIG_n
             ,data = KG4, Hess = T)

결과 요약

stargazer(reg1, reg2, reg3, type = "text")

## 
## =====================================================
##                              Dependent variable:     
##                         -----------------------------
##                         SEXATT1_n SEXATT2_n SEXATT3_n
##                            (1)       (2)       (3)   
## -----------------------------------------------------
## SEX_nFemale             0.310***  0.465***  -0.465***
##                          (0.080)   (0.095)   (0.092) 
##                                                      
## AGE_n                   0.041***    0.007   0.023*** 
##                          (0.004)   (0.004)   (0.004) 
##                                                      
## MARITAL_n기혼              -0.051     0.061     0.147  
##                          (0.115)   (0.133)   (0.125) 
##                                                      
## MARITAL_n사별, 이혼, 별거, 동거  -0.189    -0.227    -0.029  
##                          (0.166)   (0.196)   (0.189) 
##                                                      
## party_n1                0.266***    0.124   0.516*** 
##                          (0.092)   (0.112)   (0.110) 
##                                                      
## party_n2                  0.035     0.083     0.118  
##                          (0.097)   (0.115)   (0.107) 
##                                                      
## PARTYLR_n               0.156***   0.00001  0.141*** 
##                          (0.040)   (0.048)   (0.047) 
##                                                      
## EDUC_n                  -0.083*** -0.109***  -0.038  
##                          (0.031)   (0.037)   (0.037) 
##                                                      
## FAMILY_n2               0.222***  0.156***  0.435*** 
##                          (0.040)   (0.049)   (0.048) 
##                                                      
## YEAR_n2                 -1.198***   0.114   -0.843***
##                          (0.127)   (0.144)   (0.135) 
##                                                      
## RELIG_n불교                -0.106     0.003     0.102  
##                          (0.123)   (0.151)   (0.152) 
##                                                      
## RELIG_n개신교              0.438***  0.425***  0.788*** 
##                          (0.121)   (0.152)   (0.156) 
##                                                      
## RELIG_n천주교               0.358**    0.266     0.330  
##                          (0.176)   (0.228)   (0.212) 
##                                                      
## YEAR_n2:RELIG_n불교        0.361*    -0.333     0.086  
##                          (0.197)   (0.230)   (0.219) 
##                                                      
## YEAR_n2:RELIG_n개신교        0.312    -0.261    0.433*  
##                          (0.199)   (0.244)   (0.239) 
##                                                      
## YEAR_n2:RELIG_n천주교       -0.300    -0.247    -0.120  
##                          (0.268)   (0.326)   (0.290) 
##                                                      
## -----------------------------------------------------
## Observations              2,539     2,539     2,539  
## =====================================================
## Note:                     *p<0.1; **p<0.05; ***p<0.01

해석:

여성은 혼전, 혼외 성관계에 대한 부정적 태도가 남성보다 강하지만, 동성간 성관계는 남성이 더 부정적이다.
혼외 성관계를 제외한 나머지 두 유형의 행위에 대해서 연령이 높을수록 부정적이다.
국힘당 계열 지지자들은 비지지자에 비해서 혼전, 동성간 성관계에 대해서 부정적이다.
정치적으로 보수적일수록 혼전, 동성간 성관계에 대해서 부정적이다.
교육수준이 높을수록 혼전, 혼외 성관계에 대해서 부정적이지 않다.
성역할규범이 강할수록 혼전, 혼외, 동성간 성관계에 대해서 부정적이다.
개신교는 세 행위 모두에 대해서 다른 종교에 비해서 부정적이다. 또한 천주교는 혼전 성관계에 대해서 부정적이다.
년도가 지남에 따라서 혼전, 동성간 성관계에 대한 부정적 태도는 줄어든다.
불교는 년도가 지나도 혼전 성관계에 대한 부정적 태도가 줄어드는 경향이 약하다.
개신교는 년도가 지나도 동성간 성관계에 대한 부정적 태도가 줄어드는 경향이 약하다.

Step 3. 각 결과 ODDS RATIO 표

### 이거 정말 고생고생했음… Odds Ratio를 그림으로 나타낸 것인데 그림에서 노란 원의 경우 0.05 수준에서 유의한 변수들이다.

혼전 성관계

ctable1 <- coef(summary(reg1))
p <- pnorm(abs(ctable1[, "t value"]), lower.tail = FALSE) * 2
ctable1 <- cbind(ctable1, "p_value_1" = p) # 오즈 + p value
ctable2 <- exp(cbind(OR_1 = coef(reg1), confint(reg1))) #오즈비 + 신뢰구간 표

## Waiting for profiling to be done...

ctable1 <- data.frame(ctable1) #data frame으로
ctable2 <- data.frame(ctable2) #data frame으로
ctable1 <- tibble::rownames_to_column(ctable1, "boxLabels") #row이름을 1열로 + 변수명 부여
ctable2 <- tibble::rownames_to_column(ctable2, "boxLabels") #row이름을 1열로 + 변수명 부여
ctable3 <- merge(ctable1, ctable2, by = "boxLabels")
ctable3$p <- factor(ifelse(ctable3$p_value_1 < 0.05, 1, 0)) #유의하면 1 아니면 0
# ctable이라는 결과표 data frame을 구성함

ggplot(ctable3, aes(x = OR_1, y = boxLabels, shape = p)) + 
    geom_vline(aes(xintercept = 1), size = .25, linetype = "dashed") + 
    geom_errorbarh(aes(xmax = X97.5.., xmin = X2.5..), size = .5
                   , height = .2, color = "gray50") +
    geom_point(col = "orange", size = 3.0) +
    scale_shape_manual(values=c(1, 16))+
    theme_bw()+
    theme(panel.grid.minor = element_blank()) +
    ylab("") +
    xlab("Odds ratio") +
    ggtitle("") +
    scale_x_continuous(breaks = 0.5*(1:24)) +
    theme(legend.position = "none")

혼전 성관계 해석:

1순위: 년도가 지남에 따라서 부정적 태도는 약 1.7배 낮아진다

2순위: 개신교는 무교를 기준으로 볼때 약 1.5배 부정적이다 기타 등등…

혼외 성관계

ctable1 <- coef(summary(reg2))
p <- pnorm(abs(ctable1[, "t value"]), lower.tail = FALSE) * 2
ctable1 <- cbind(ctable1, "p_value_2" = p) # 오즈 + p value
ctable2 <- exp(cbind(OR_2 = coef(reg2), confint(reg2))) #오즈비 + 신뢰구간 표

## Waiting for profiling to be done...

ctable1 <- data.frame(ctable1) #data frame으로
ctable2 <- data.frame(ctable2) #data frame으로
ctable1 <- tibble::rownames_to_column(ctable1, "boxLabels") #row이름을 1열로 + 변수명 부여
ctable2 <- tibble::rownames_to_column(ctable2, "boxLabels") #row이름을 1열로 + 변수명 부여
ctable4 <- merge(ctable1, ctable2, by = "boxLabels")
ctable4$p <- factor(ifelse(ctable4$p_value_2 < 0.05, 1, 0)) #유의하면 1 아니면 0
# ctable이라는 결과표 data frame을 구성함

ggplot(ctable4, aes(x = OR_2, y = boxLabels, shape = p)) + 
    geom_vline(aes(xintercept = 1), size = .25, linetype = "dashed") + 
    geom_errorbarh(aes(xmax = X97.5.., xmin = X2.5..), size = .5
                   , height = .2, color = "gray50") +
    geom_point(col = "orange", size = 3.0) +
    scale_shape_manual(values=c(1, 16))+
    theme_bw()+
    theme(panel.grid.minor = element_blank()) +
    ylab("") +
    xlab("Odds ratio") +
    ggtitle("") +
    scale_x_continuous(breaks = 0.5*(1:24)) +
    theme(legend.position = "none")

혼외 성관계 해석:

혼외 성관계의 경우 년도가 바뀌어도 대다수가 절대적으로 부정적이다.

따라서 혼전 성관계와 비교할 때 유의한 요인들의 수가 적다.

특히 흥미로운 것은 여성이 남성에 비해서 약 1.6배정도 더 혼외 성관계에 부정적이라는 점이다.

기타 등등…

동성간 성관계

ctable1 <- coef(summary(reg3))
p <- pnorm(abs(ctable1[, "t value"]), lower.tail = FALSE) * 2
ctable1 <- cbind(ctable1, "p_value_3" = p) # 오즈 + p value
ctable2 <- exp(cbind(OR_3 = coef(reg3), confint(reg3))) #오즈비 + 신뢰구간 표

## Waiting for profiling to be done...

ctable1 <- data.frame(ctable1) #data frame으로
ctable2 <- data.frame(ctable2) #data frame으로
ctable1 <- tibble::rownames_to_column(ctable1, "boxLabels") #row이름을 1열로 + 변수명 부여
ctable2 <- tibble::rownames_to_column(ctable2, "boxLabels") #row이름을 1열로 + 변수명 부여
ctable5 <- merge(ctable1, ctable2, by = "boxLabels")
ctable5$p <- factor(ifelse(ctable5$p_value_3 < 0.05, 1, 0)) #유의하면 1 아니면 0
# ctable이라는 결과표 data frame을 구성함

ggplot(ctable5, aes(x = OR_3, y = boxLabels, shape = p)) + 
    geom_vline(aes(xintercept = 1), size = .25, linetype = "dashed") + 
    geom_errorbarh(aes(xmax = X97.5.., xmin = X2.5..), size = .5
                   , height = .2, color = "gray50") +
    geom_point(col = "orange", size = 3.0) +
    scale_shape_manual(values=c(1, 16))+
    theme_bw()+
    theme(panel.grid.minor = element_blank()) +
    ylab("") +
    xlab("Odds ratio") +
    ggtitle("") +
    scale_x_continuous(breaks = 0.5*(1:24)) +
    theme(legend.position = "none")

동성간 성관계 해석:

동성간 성관계에 대한 태도에서는 역시나 개신교의 효과가 크게 나타난다.

그럼에도 년도가 지남에 따라서 부정적 태도는 감소하고 있다.

앞선 두 행위와 달리 여성은 남성에 비해서 약 1.4배 우호적인 태도를 보인다. 기타 등등…

Step 4. 종교의 효과(Marginal Effect at Mean)

혼전 성관계

model1 <- ggpredict(reg1, terms = c("RELIG_n"))
ggplot(model1) +
    aes(x = x, fill = x, weight = predicted) +
    geom_bar(width = 0.5) +
    geom_point(aes(y = predicted), size = 1.5)+
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), col = 1, alpha = 0.5, width = 0.3)+
    scale_fill_hue() +
    theme_bw()+
    labs(title = "혼전 성관계 부정적 태도"
         , subtitle = ""
         , x = "RELIGION"
         , y = "Predicted"
         , fill = "RELIGION")  + 
    facet_grid(. ~ response.level)

혼외 성관계

model2 <- ggpredict(reg2, terms = c("RELIG_n"))
ggplot(model2) +
    aes(x = x, fill = x, weight = predicted) +
    geom_bar(width = 0.5) +
    geom_point(aes(y = predicted), size = 1.5)+
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), col = 1, alpha = 0.5, width = 0.3)+
    scale_fill_hue() +
    theme_bw()+
    labs(title = "혼외 성관계 부정적 태도"
         , subtitle = ""
         , x = "RELIGION"
         , y = "Predicted"
         , fill = "RELIGION")  + 
    facet_grid(. ~ response.level)

동성간 성관계

model3 <- ggpredict(reg3, terms = c("RELIG_n"))
ggplot(model3) +
    aes(x = x, fill = x, weight = predicted) +
    geom_bar(width = 0.5) +
    geom_point(aes(y = predicted), size = 1.5)+
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), col = 1, alpha = 0.5, width = 0.3)+
    scale_fill_hue() +
    theme_bw()+
    labs(title = "동성간 성관계 부정적 태도"
         , subtitle = ""
         , x = "RELIGION"
         , y = "Predicted"
         , fill = "RELIGION")  + 
    facet_grid(. ~ response.level)

Step 5. 종교 x 년도의 효과

혼전 성관계

model4 <- ggpredict(reg1, terms = c("YEAR_n2", "RELIG_n"))
ggplot(model4) +
    aes(x = x, y = predicted, colour = group) +
    geom_line(size = 1L) +
    scale_color_hue() +
    theme_bw()+ 
    scale_x_continuous(breaks = c(0, 1), labels=c("08년", "18년"))+
    labs(title = "혼전 성관계 부정적 태도 : 년도*종교"
         , subtitle = ""
         , x = "YEAR"
         , y = "Probability"
         , col = "Religion") +
facet_grid(. ~ response.level)

혼외 성관계

model5 <- ggpredict(reg2, terms = c("YEAR_n2", "RELIG_n"))
ggplot(model5) +
    aes(x = x, y = predicted, colour = group) +
    geom_line(size = 1L) +
    scale_color_hue() +
    theme_bw()+ 
    scale_x_continuous(breaks = c(0, 1), labels=c("08년", "18년"))+
    labs(title = "혼외 성관계 부정적 태도 - 년도*종교"
         , subtitle = ""
         , x = "YEAR"
         , y = "Probability"
         , col = "Religion") +
    facet_grid(. ~ response.level)

동성간 성관계

model6 <- ggpredict(reg3, terms = c("YEAR_n2", "RELIG_n"))
ggplot(model6) +
    aes(x = x, y = predicted, colour = group) +
    geom_line(size = 1L) +
    scale_color_hue() +
    theme_bw()+ 
    scale_x_continuous(breaks = c(0, 1), labels=c("08년", "18년"))+
    labs(title = "동성 성관계 부정적 태도 - 년도*종교"
         , subtitle = ""
         , x = "YEAR"
         , y = "Probability"
         , col = "Religion") +
    facet_grid(. ~ response.level)

동성간 성관계 中 Tatally Wrong만 보기

# 동성 성관계 종교x년도 

model7 <- ggpredict(reg3, terms = c("YEAR_n2", "RELIG_n"))
model7 %>% filter(response.level == c("4_Totally Wrong")) %>%
ggplot() +
    aes(x = x, y = predicted, colour = group) +
    geom_line(size = 1L) +
    scale_color_hue() +
    theme_bw()+ 
    scale_x_continuous(breaks = c(0, 1), labels=c("08년", "18년"))+
    labs(title = "동성 성관계 부정적 태도 - 년도*종교"
         , subtitle = ""
         , x = "YEAR"
         , y = "Probability"
         , col = "Religion") +
    facet_grid(. ~ response.level)

한 눈에 보기에도 가장 두드러지는 것이 개신교인들의 동성애에 대한 태도 변화다. 물론 년도가 지남에 따라서 개신교인들 역시 부정적 태도가 줄어들고 있긴 하다. 하지만 그 정도가 다른 종교에 비해서 확연히 완만하다.

따라서 앞으로의 작업은 왜 개신교인들만 이렇게 변화가 느릴까?에 대한 답을 찾는 것이다.

사실 지금 작업중이긴 한데…그건 논문으로 써야해서!

아 키보드 좋다!

‘비규범적’ 성관계에 대한 종교인의 태도

HeeYoung

2021 3 14

데이터

Step 1. 변수

종속변수는 각각 혼전 성관계, 혼외 성관계, 동성간 성관계에 대한 4점짜리 척도다. 해당 행위가 얼마나 “잘못되었는가”에 대한 답변이다.

* 혼전 성관계 각각 4첨 척도 빈도

* 혼외 성관계 각각 4첨 척도 빈도

* 동성간 성관계 각각 4첨 척도 빈도

설명변수는 1)성별(SEX_n) 2)연령(AGE_n) 3)교육수준(EDUC_n) 4)정치성향(PARTYLR_n) 5)혼인(MARITAL_n) 6)종교(RELIG_n) 6)년도(YEAR_n2) 7)국힘당 지지(party_n1) 8)더민주 지지(party_n2) 9)성역할규범(FAMILY_n2)이다.

혼외 성관계

동성간 성관계

결과 요약

해석:

Step 3. 각 결과 ODDS RATIO 표

### 이거 정말 고생고생했음… Odds Ratio를 그림으로 나타낸 것인데 그림에서 노란 원의 경우 0.05 수준에서 유의한 변수들이다.

혼전 성관계

혼외 성관계

동성간 성관계

Step 4. 종교의 효과(Marginal Effect at Mean)

혼전 성관계

혼외 성관계

동성간 성관계

Step 5. 종교 x 년도의 효과

혼전 성관계

혼외 성관계

동성간 성관계

동성간 성관계 中 Tatally Wrong만 보기