통계상담과제

통계상담과제는 실제로 학생들에게 설문조사를 하여 분석한 결과입니다.

-------------------------------------------------------------------------------------------------------------------------------

getwd()
setwd("G:\\통계상담")
setwd("/media/sukjune/SANDISK/통계상담")

save.image("아침식사.RData") # 변수저장
load("아침식사.RData") # 변수 불러오기

breakfast <- read.csv(file = "아침식사자료(수정).csv", stringsAsFactors = FALSE, fileEncoding = "UTF-8", colClasses = "character")
head(breakfast) # 아침식사데이터

# ------------------------------------------- 1. 데이터 변형 --------------------------------------------------------
# install.packages("dplyr")
library(dplyr)

str(breakfast)
names(breakfast) <- tolower(names(breakfast))
breakfast <- breakfast %>% mutate(
x2 = as.numeric(x2), # 키 숫자로 변환
x3 = as.numeric(x3) # 몸무게 숫자로 변환
) %>% select(-c(x, time, x16, x17)) # 관련없는 변수 제거

summary(breakfast) # 요약 통계량
# x3에 결측존재
# x12에 응답안한 사람 ERROR로 표시
# x18에 공백존재

delete_list <- complete.cases(breakfast) # 몸무게의 결측 위치
breakfast <- breakfast[delete_list,] # 결측제거

summary(breakfast)
# x3의 결측제거
# x18의 공백 제거 - 추후 level을 변경해야함
# "#N/A"를 "해당없음"로 바꿔야한다.

breakfast[breakfast$x12 == "#N/A","x12"] <- "해당없음" # "#N/A"를 "해당없음"으로 바꿈

x <- breakfast %>% select(-x15) # 설명변수들의 데이터프레임
names(x) <- paste("x",1:18, sep = "") # 이름 재정의
y <- breakfast$x15

rm(list = setdiff(ls(), c("breakfast","y","x")))
# ---------------------------------- 2. 각 변수별 관계 --------------------------------------------------
# 유의수준 0.1로 설정
# ---- 1) 성별과 아침식사 ----
a <- table(x$x1, y, dnn = list("성별", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x1, y) # 독립성 검정
# p-value가 유의수준 0.1보다 작지 않으므로 서로 독립이다.

# ---- 2) 학년과 아침식사(x) ----
a <- table(x$x4, y, dnn = list("학년", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x4, y) # 독립성 검정
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 3) 흡연과 아침식사 ----
a <- table(x$x5, y, dnn = list("흡연", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x5, y) # 독립성 검정
# p-value가 유의수준 0.1보다 크므로 서로 독립이다.

# ---- 4) 수면시간과 아침식사 ----
table(x$x6, y, dnn = list("수면시간", "아침식사"))
# 3-4시간의 관측치가 작아 "4-5시간","5-6시간"과 통합한다.

x2 <- x %>% mutate(
x6 = ifelse(x6 %in% c("3-4시간 미만", "4-5시간 미만", "5-6시간 미만"),"6시간 미만",x6) # 통합
)
# "7시간 이상"으로 통합한다.
x2 <- x2 %>% mutate(
x6 = ifelse(x6 %in% c("7-8시간 미만","8시간 이상"),"7시간 이상",x6) # 통합
)
a <- table(x2$x6, y, dnn = list("수면시간", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x6, y) # 범주에서 관측치가 적어 통합하여 분석
# p-value가 유의수준 0.1보다 크므로 서로 독립이다.

# ---- 5) 음주과 아침식사 ----
table(x$x7, y, dnn = list("음주", "아침식사"))
# 주 4회이상으로 통합한다.

x2$x7 <- ifelse(x$x7 %in% c("주 4-5회", "주 6-7회"), "주 4회이상",x$x7)
a <- table(x2$x7, y, dnn = list("음주", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x7,y)
# p-value가 유의수준 0.1보다 크므로 서로 독립이다.

# ---- 6) 운동과 아침식사(x) ----
table(x$x8, y, dnn = list("운동", "아침식사"))
# "주 4회이상"으로 통합한다.

x2$x8 <- ifelse(x$x8 %in% c("주 4-5회", "주 6-7회"),"주 4회이상",x$x8)
a <- table(x2$x8, y, dnn = list("운동", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x8, y)
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 7) 야식과 아침식사(x) ----
table(x$x9, y, dnn = list("야식", "아침식사"))
# "주 4회이상"으로 합친다.

x2$x9 <- ifelse(x$x9 %in% c("주 4-5회", "주 6-7회"),"주 4회이상",x$x9)
a <- table(x2$x9, y, dnn = list("야식", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x9, y) # 정확하지 않을 수 있어 fisher.test를사용
fisher.test(x2$x9, y, alternative = "two.sided")
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 8) 간식과 아침식사 ----
table(x$x10, y, dnn = list("간식", "아침식사"))
# "주 4회이상"으로 합한다.

x2$x10 <- ifelse(x$x10 %in% c("주 4-5회", "주 6-7회"),"주 4회이상",x$x10)
a <- table(x2$x10, y, dnn = list("간식", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x10, y)
# p-value가 유의수준 0.1보다 크므로 서로 독립이다.

# ---- 9) 주거형태과 아침식사(x) ----
a <- table(x$x11, y, dnn = list("주거형태", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x11, y)
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 10) 소요시간과 아침식사(x) ----
table(x$x12, y, dnn = list("소요시간", "아침식사"))
# "30분 미만"과 "30분 이상"으로 구분

x2$x12 <- ifelse(x$x12 == "해당없음", "해당없음",
ifelse(x$x12 %in% c("10분 미만", "10분 이상 30분 미만"),"30분 미만", "30분 이상"))
a <- table(x2$x12, y, dnn = list("소요시간", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x12, y) # 정확검정을 시행
fisher.test(x2$x12, y, alternative = "two.sided")
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 11) 수업시간과 아침식사(x) ----
table(x$x13, y, dnn = list("수업시간", "아침식사"))
# "12시간 이상 18시간 미만"으로 통합

x2$x13 <- ifelse(x$x13 %in% c("12시간 이상 15시간 미만", "15시간 이상 18시간 미만"),
"12시간 이상 18시간 미만", x$x13)
a <- table(x2$x13, y, dnn = list("수업시간", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x13, y) # 정확검정 실시
fisher.test(x2$x13, y, alternative = 'two.sided')
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 12) 등교횟수과 아침식사(x) ----
table(x$x14, y, dnn = list("등교횟수", "아침식사"))
# "주 3회이하"로 통합한다.

x2$x14 <- ifelse(x$x14 %in% c("주 2회", "주 3회"),"주 3회이하",x$x14)
a <- table(x2$x14, y, dnn = list("등교횟수", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x2$x14, y)
fisher.test(x2$x14, y, alternative = 'two.sided')
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 13) 아침-건강는과 아침식사(x) ----
# 아침식사가 건강에 영향을 미치는지에 대한 문항
table(x$x15, y, dnn = list("아침-건강", "아침식사"))
# "별 관련없다"로 통합

x2$x15 <- ifelse(x$x15 %in% c("관련은 있으나 큰 영향을 미치지는 않는다","아침식사와 건강과는 별 관련이 없다"),
"관련없다.",x$x15)
a <- table(x2$x15, y, dnn = list("아침-건강", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x15, y) # 정확검정 실시
fisher.test(x$x15, y, alternative = "two.sided", workspace = 1e+8)
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 14) 아침인식과 아침식사(x) ----
a <- table(x$x16, y, dnn = list("아침인식", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x16, y)
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 15) 점심인식과 아침식사(x) ----
a <- table(x$x17, y, dnn = list("점심인식", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x17, y)
fisher.test(x$x17, y, alternative = "two.sided")
# p-value가 유의수준 0.1보다 작으므로 서로 독립이 아니다.

# ---- 16) 저녁인식과 아침식사 ----
a <- table(x$x18, y, dnn = list("저녁인식", "아침식사"))
addmargins(a, margin =2)
round(prop.table(a,margin = 1),2)
chisq.test(x$x18, y)
# p-value가 유의수준 0.1보다 크므로 서로 독립이다.

x <- x2 %>% select(x4, x8, x9, x11, x12, x13, x14, x15, x16, x17)
rm(list = setdiff(ls(), c("breakfast", "x", "y", "x2")))
# 최종적으로 x4, x8, x9, x11, x12, x13, x14, x15, x16, x17의 변수만 이용하여 로지스틱 회귀분석을 실시한다.

# ---------------------------------- 3. 로짓스틱 회귀모형 --------------------------------------------------
# 널 모형에서 모든 교호작용이 존재하는 모형까지 전진선택법으로 변수를 선택한다.
# 유의수준은 0.1로 사용
# 모형은 다범주 로지스틱을 이용한다.

# ---- 1) 데이터 전처리 ----
# x11과 x12를 통합하여 s를 만든다.
x$s <- ifelse(x$x11 == "통학",x$x12,x$x11)

x <- x %>% mutate(
x4 = factor(x4, levels = c("1학년", "2학년", "3학년", "4학년")),
x8 = factor(x8, levels = c("주 0-1회", "주 2-3회", "주 4회이상")),
x9 = factor(x9, levels = c("주 0-1회", "주 2-3회", "주 4회이상")),
x11 = factor(x11, levels = c("기숙사", "자취", "통학")),
x12 = factor(x12, levels = c("30분 미만", "30분 이상", "해당없음")),
x13 = factor(x13, levels = c("12시간 미만", "12시간 이상 18시간 미만", "18시간 이상")),
x14 = factor(x14, levels = c("주 3회이하", "주 4회", "주 5회")),
x15 = factor(x15, levels = c("잘 모르겠다", "관련없다.", "아침식사는 건강에 큰 영향을 미친다")),
x16 = factor(x16, levels = c("걸러도 된다", "보통", "먹어야 한다")),
x17 = factor(x17, levels = c("걸러도 된다", "보통", "먹어야 한다")),
s = factor(s, levels = c("기숙사", "자취", "30분 미만", "30분 이상"))
)
str(x)
summary(x)

# ---- 2)변수선택 ----
# install.packages("nnet")
library(nnet)
y <- factor(y)
model_null <- multinom(y ~ 1, data = x, family = binomial("logit")) # 널 모형

result <- step(model_null, scope = list(lower = ~ 1, upper = ~ x4*x8*x9*s*x13*x14*x15*x16*x17), direction = "forward")
var_f <- formula(result) # 선택된 함수식

# ---- 3) 모형생성 ----
model <- multinom(var_f, data = x, family = binomial("logit"))

# ---- 4) 유의성 검정 ----
anova(model_null, model) # 모형의 유의성검정

# 회귀계수의 유의성 검정
summary(model) # p-value를 직접 계산해야한다.
z_static <- summary(model)$coefficients/summary(model)$standard.errors # 표준 정규분포를 따른다.
p_value <- (1 -pnorm(abs(z_static), mean = 0, sd = 1))*2
p_logical <- p_value > 0.1; p_logical # 유의수준 0.1보다 크면 TRUE
# x16 - "보통"과 "걸러도 된다" -> "걸러도 된다"로 통합
# s - "자취"와 "기숙사"를 통합
# x15제외

# ---- 5) 데이터 전처리2 ----
x <- x2 %>% select(x4, x8, x9, x11, x12, x13, x14, x15, x16, x17)
x$s <- ifelse(x$x11 == "통학",x$x12,"기숙,자취")
x$x16 <- ifelse(x$x16 %in% c("보통", "걸러도 된다"),"걸러도 된다","먹어야 한다")

x <- x %>% mutate(
x4 = factor(x4, levels = c("1학년", "2학년", "3학년", "4학년")),
x8 = factor(x8, levels = c("주 0-1회", "주 2-3회", "주 4회이상")),
x9 = factor(x9, levels = c("주 0-1회", "주 2-3회", "주 4회이상")),
x11 = factor(x11, levels = c("기숙사", "자취", "통학")),
x12 = factor(x12, levels = c("30분 미만", "30분 이상", "해당없음")),
x13 = factor(x13, levels = c("12시간 미만", "12시간 이상 18시간 미만", "18시간 이상")),
x14 = factor(x14, levels = c("주 3회이하", "주 4회", "주 5회")),
x15 = factor(x15, levels = c("잘 모르겠다", "관련없다.", "아침식사는 건강에 큰 영향을 미친다")),
x16 = factor(x16, levels = c("걸러도 된다", "먹어야 한다")),
x17 = factor(x17, levels = c("걸러도 된다", "보통", "먹어야 한다")),
s = factor(s, levels = c("기숙,자취", "30분 미만", "30분 이상"))
)
str(x)
summary(x)

# ---- 6) 모형생성2 ----
model <- multinom(var_f, data = x)
model <- update(model, ~ . -x15) # x15제거, 절편제거

# ---- 7) 유의성 검정2 ----
anova(model_null, model) # 모형의 유의성검정

# 회귀계수의 유의성 검정
summary(model) # p-value를 직접 계산해야한다.
z_static <- summary(model)$coefficients/summary(model)$standard.errors # 표준 정규분포를 따른다.
p_value <- (1 -pnorm(abs(z_static), mean = 0, sd = 1))*2
p_logical <- p_value > 0.1; p_logical # 유의수준 0.1보다 크면 TRUE

# ---- 오즈 ----
# 오즈는 0-1회일 확률과의 비율이다.
Beta <- coef(model)
Beta[p_logical] <- 0 # 유의하지 않는 회귀계수는 0으로 바꾼다.
Beta
OR <- exp(Beta)[,-1]
OR
# "1학년"
# x16 - "먹어야 한다"라고 생각해도 실제로 먹지 않는다는 것을 알 수 있다.
# s - "통학"하는 사람은 자취나 기숙사에 사는 사람보다 아침을 먹지는것을 알 수 있다.
# x14 - 자주 학교에 등교하는 사람들이 아침을 잘 먹는것을 알 수 있다.

Beta2 <- Beta[-1,]
Beta2 <- Beta2 - rbind(Beta[1,],Beta[1,]) # "주 2-3회"가 분모인 회귀계수
OR2 <- exp(Beta2)[,-1]
OR2

Beta3 <- Beta2[2,] - Beta2[1,] # "주 4-5회"가 분모인 회귀계수
OR3 <- exp(Beta3)[-1]
OR3