R로 EFA하기
1. 자료 불러오기
setwd("C:/Users/LG/Desktop/efa")
#setwd = R에서 자료를 불러오려면 자료가 어디에 저장되어있는지 설정해야함.
#setwd로 지정한 폴더로 R 자료가 저장되고 폴더내에 있는 파일을 불러올 수 있음.
#나는 바탕화면에 자료를 저장해놔서 저장루트가 위와 같음.
library(readxl)
## Warning: package 'readxl' was built under R version 3.5.3
#xlsx 확장자를 가지고 있는 엑셀파일을 불러오려면 readxl이라는 패키지를 불러와서 실행해야함.
data <- read_excel(path="depression.xlsx", sheet="depress", col_names=TRUE)
#우울척도 파일을 불러와서 'data'로 저장함.
head(data)
## # A tibble: 6 x 15
## B6_1 B6_2 B6_3 B6_4 B6_5 B6_6 B6_7 B6_8 B6_9 B6_10 B6_11 B6_12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 1 2 1 2 2 2 1 2 1 2 1
## 2 1 1 2 2 1 2 1 2 2 2 1 1
## 3 1 1 2 2 1 2 1 2 2 2 1 1
## 4 2 2 1 2 2 1 2 2 2 2 2 2
## 5 1 1 2 1 1 2 1 2 2 2 1 2
## 6 1 1 2 2 1 2 1 2 2 1 1 2
## # ... with 3 more variables: B6_13 <dbl>, B6_14 <dbl>, B6_15 <dbl>
#'data'로 저장된 우울척도 파일을 처음부터 5개만 보여줌.
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3
##
## 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
# %>% 연산자를 사용하기 위해서 불러온 패키지
data %>% dim()
## [1] 10299 15
data %>% str()
## Classes 'tbl_df', 'tbl' and 'data.frame': 10299 obs. of 15 variables:
## $ B6_1 : num 2 1 1 2 1 1 2 2 2 2 ...
## $ B6_2 : num 1 1 1 2 1 1 1 1 1 1 ...
## $ B6_3 : num 2 2 2 1 2 2 2 2 2 2 ...
## $ B6_4 : num 1 2 2 2 1 2 1 1 2 1 ...
## $ B6_5 : num 2 1 1 2 1 1 1 1 2 2 ...
## $ B6_6 : num 2 2 2 1 2 2 2 1 1 1 ...
## $ B6_7 : num 2 1 1 2 1 1 1 1 2 2 ...
## $ B6_8 : num 1 2 2 2 2 2 2 2 2 2 ...
## $ B6_9 : num 2 2 2 2 2 2 2 1 2 1 ...
## $ B6_10: num 1 2 2 2 2 1 1 1 1 2 ...
## $ B6_11: num 2 1 1 2 1 1 1 1 1 2 ...
## $ B6_12: num 1 1 1 2 2 2 2 2 2 2 ...
## $ B6_13: num 2 1 1 2 2 1 2 1 2 2 ...
## $ B6_14: num 2 2 2 2 2 2 2 2 2 2 ...
## $ B6_15: num 2 2 2 2 2 2 2 2 2 2 ...
data %>% summary()
## B6_1 B6_2 B6_3 B6_4
## Min. :1.00 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.00 1st Qu.:1.000 1st Qu.:2.000 1st Qu.:1.000
## Median :1.00 Median :2.000 Median :2.000 Median :2.000
## Mean :1.53 Mean :1.737 Mean :1.992 Mean :1.831
## 3rd Qu.:2.00 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :9.00 Max. :9.000 Max. :9.000 Max. :9.000
## B6_5 B6_6 B6_7 B6_8
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:2.000 1st Qu.:1.000 1st Qu.:2.000
## Median :1.000 Median :2.000 Median :1.000 Median :2.000
## Mean :1.577 Mean :1.988 Mean :1.539 Mean :2.006
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :9.000 Max. :9.000 Max. :9.000 Max. :9.000
## B6_9 B6_10 B6_11 B6_12
## Min. :1.000 Min. :1.00 Min. :1.000 Min. :1.000
## 1st Qu.:2.000 1st Qu.:1.00 1st Qu.:1.000 1st Qu.:2.000
## Median :2.000 Median :2.00 Median :1.000 Median :2.000
## Mean :1.998 Mean :1.88 Mean :1.514 Mean :1.965
## 3rd Qu.:2.000 3rd Qu.:2.00 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :9.000 Max. :9.00 Max. :9.000 Max. :9.000
## B6_13 B6_14 B6_15
## Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:2.000 1st Qu.:2.000
## Median :1.000 Median :2.000 Median :2.000
## Mean :1.616 Mean :1.959 Mean :1.992
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :9.000 Max. :9.000 Max. :9.000
#data 전체 n 수와 변수 개수를 파악하는 명령어 dim
#data 속성을 확인하기 위한 명령어 str
#data 전체 요약을 위해서 사용하는 명령어 summary
2. 데이터 수정하기
data[data == 9]<-NA
#summary에서 9가 있으니 9는 NA로 변환
cdata <- na.omit(data)
#NA값을 가지는 자료는 모두 지우기
head(cdata)
## # A tibble: 6 x 15
## B6_1 B6_2 B6_3 B6_4 B6_5 B6_6 B6_7 B6_8 B6_9 B6_10 B6_11 B6_12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 1 2 1 2 2 2 1 2 1 2 1
## 2 1 1 2 2 1 2 1 2 2 2 1 1
## 3 1 1 2 2 1 2 1 2 2 2 1 1
## 4 2 2 1 2 2 1 2 2 2 2 2 2
## 5 1 1 2 1 1 2 1 2 2 2 1 2
## 6 1 1 2 2 1 2 1 2 2 1 1 2
## # ... with 3 more variables: B6_13 <dbl>, B6_14 <dbl>, B6_15 <dbl>
cdata[cdata==1]<-1 #예
cdata[cdata==2]<-0 #아니오
#1,2 값을 1,0으로 바꾸기
rcdata <- cdata %>% mutate(reB6_1 = ifelse(B6_1 == 1,0,1),
reB6_5 = ifelse(B6_5 == 1,0,1),
reB6_7 = ifelse(B6_7 == 1,0,1),
reB6_11 = ifelse(B6_11 == 1,0,1),
reB6_13 = ifelse(B6_13 == 1,0,1))
#역문항변환
ddd <- cdata[c(1:5000),]
write.csv(ddd, "ddd.csv")
data1 <- rcdata[c(1:5000),]
data2 <- rcdata[c(5001:10000),]
#자료나누기
head(data1)
## # A tibble: 6 x 20
## B6_1 B6_2 B6_3 B6_4 B6_5 B6_6 B6_7 B6_8 B6_9 B6_10 B6_11 B6_12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 1 0 1 0 0 0 1 0 1 0 1
## 2 1 1 0 0 1 0 1 0 0 0 1 1
## 3 1 1 0 0 1 0 1 0 0 0 1 1
## 4 0 0 1 0 0 1 0 0 0 0 0 0
## 5 1 1 0 1 1 0 1 0 0 0 1 0
## 6 1 1 0 0 1 0 1 0 0 1 1 0
## # ... with 8 more variables: B6_13 <dbl>, B6_14 <dbl>, B6_15 <dbl>,
## # reB6_1 <dbl>, reB6_5 <dbl>, reB6_7 <dbl>, reB6_11 <dbl>, reB6_13 <dbl>
head(data2)
## # A tibble: 6 x 20
## B6_1 B6_2 B6_3 B6_4 B6_5 B6_6 B6_7 B6_8 B6_9 B6_10 B6_11 B6_12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 0 0 1 0 0 0 0 0 1 0
## 2 1 1 0 0 0 0 0 0 0 0 0 0
## 3 1 1 0 1 1 0 1 0 0 0 1 0
## 4 1 0 0 0 1 0 1 0 0 0 1 0
## 5 1 1 0 0 1 0 1 0 0 0 1 0
## 6 0 1 1 1 0 1 0 1 1 1 0 1
## # ... with 8 more variables: B6_13 <dbl>, B6_14 <dbl>, B6_15 <dbl>,
## # reB6_1 <dbl>, reB6_5 <dbl>, reB6_7 <dbl>, reB6_11 <dbl>, reB6_13 <dbl>
data1 <- data1 %>% select(-c(B6_1,B6_5,B6_7,B6_11,B6_13))
data2 <- data2 %>% select(-c(B6_1,B6_5,B6_7,B6_11,B6_13))
#자료나누기 후에 역문항에 해당하는 문항 제외
head(data1)
## # A tibble: 6 x 15
## B6_2 B6_3 B6_4 B6_6 B6_8 B6_9 B6_10 B6_12 B6_14 B6_15 reB6_1 reB6_5
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 1 0 1 0 1 1 0 0 1 1
## 2 1 0 0 0 0 0 0 1 0 0 0 0
## 3 1 0 0 0 0 0 0 1 0 0 0 0
## 4 0 1 0 1 0 0 0 0 0 0 1 1
## 5 1 0 1 0 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 0 1 0 0 0 0 0
## # ... with 3 more variables: reB6_7 <dbl>, reB6_11 <dbl>, reB6_13 <dbl>
head(data2)
## # A tibble: 6 x 15
## B6_2 B6_3 B6_4 B6_6 B6_8 B6_9 B6_10 B6_12 B6_14 B6_15 reB6_1 reB6_5
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 0 0 0 0 1 1 0 0
## 2 1 0 0 0 0 0 0 0 1 1 0 1
## 3 1 0 1 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0
## 5 1 0 0 0 0 0 0 0 1 0 0 0
## 6 1 1 1 1 1 1 1 1 1 1 1 1
## # ... with 3 more variables: reB6_7 <dbl>, reB6_11 <dbl>, reB6_13 <dbl>
3. EFA
library(psych)
#EFA하기 위한 패키지
#1에서 5000명 데이터로 상관행렬 만들기
fdata1 <- data1 %>% mutate_all(as.numeric)
tc1 <- tetrachoric(fdata1)
tcc1 <-as.data.frame(tc1$rho)
tcc1
## B6_2 B6_3 B6_4 B6_6 B6_8 B6_9
## B6_2 1.0000000 0.6000204 0.6717808 0.5820488 0.5971534 0.5763777
## B6_3 0.6000204 1.0000000 0.6818332 0.5772471 0.7120961 0.5297541
## B6_4 0.6717808 0.6818332 1.0000000 0.5890384 0.6678793 0.6241787
## B6_6 0.5820488 0.5772471 0.5890384 1.0000000 0.7391000 0.5281833
## B6_8 0.5971534 0.7120961 0.6678793 0.7391000 1.0000000 0.5852709
## B6_9 0.5763777 0.5297541 0.6241787 0.5281833 0.5852709 1.0000000
## B6_10 0.5045364 0.4250930 0.4500735 0.4304241 0.4867057 0.4039012
## B6_12 0.5635080 0.6923963 0.6666805 0.5996274 0.6935987 0.5414342
## B6_14 0.5503661 0.7122749 0.6096022 0.6225937 0.7508038 0.5145593
## B6_15 0.5190842 0.6657157 0.5948551 0.6113043 0.7392025 0.5418847
## reB6_1 0.6169322 0.6167177 0.6122111 0.5360571 0.6205609 0.5085199
## reB6_5 0.5648171 0.5460968 0.6317238 0.5525071 0.6250587 0.5248282
## reB6_7 0.5606241 0.5421169 0.6192816 0.6027516 0.7316499 0.5332303
## reB6_11 0.4959527 0.5183079 0.5792680 0.4911839 0.6189533 0.5342828
## reB6_13 0.5644804 0.4129777 0.4994920 0.4081680 0.4771173 0.4858979
## B6_10 B6_12 B6_14 B6_15 reB6_1 reB6_5
## B6_2 0.5045364 0.5635080 0.5503661 0.5190842 0.6169322 0.5648171
## B6_3 0.4250930 0.6923963 0.7122749 0.6657157 0.6167177 0.5460968
## B6_4 0.4500735 0.6666805 0.6096022 0.5948551 0.6122111 0.6317238
## B6_6 0.4304241 0.5996274 0.6225937 0.6113043 0.5360571 0.5525071
## B6_8 0.4867057 0.6935987 0.7508038 0.7392025 0.6205609 0.6250587
## B6_9 0.4039012 0.5414342 0.5145593 0.5418847 0.5085199 0.5248282
## B6_10 1.0000000 0.5118621 0.4687222 0.4270685 0.3755534 0.3700102
## B6_12 0.5118621 1.0000000 0.7895336 0.6931435 0.5303439 0.5316572
## B6_14 0.4687222 0.7895336 1.0000000 0.8335893 0.5486260 0.5329383
## B6_15 0.4270685 0.6931435 0.8335893 1.0000000 0.6002454 0.5952835
## reB6_1 0.3755534 0.5303439 0.5486260 0.6002454 1.0000000 0.7095953
## reB6_5 0.3700102 0.5316572 0.5329383 0.5952835 0.7095953 1.0000000
## reB6_7 0.3590840 0.5667831 0.5614758 0.6119466 0.7415162 0.8730905
## reB6_11 0.4314786 0.6352459 0.5769019 0.5957207 0.6728306 0.7404991
## reB6_13 0.4950856 0.5294182 0.4922610 0.4810883 0.5098527 0.5979681
## reB6_7 reB6_11 reB6_13
## B6_2 0.5606241 0.4959527 0.5644804
## B6_3 0.5421169 0.5183079 0.4129777
## B6_4 0.6192816 0.5792680 0.4994920
## B6_6 0.6027516 0.4911839 0.4081680
## B6_8 0.7316499 0.6189533 0.4771173
## B6_9 0.5332303 0.5342828 0.4858979
## B6_10 0.3590840 0.4314786 0.4950856
## B6_12 0.5667831 0.6352459 0.5294182
## B6_14 0.5614758 0.5769019 0.4922610
## B6_15 0.6119466 0.5957207 0.4810883
## reB6_1 0.7415162 0.6728306 0.5098527
## reB6_5 0.8730905 0.7404991 0.5979681
## reB6_7 1.0000000 0.7947694 0.5988176
## reB6_11 0.7947694 1.0000000 0.5786799
## reB6_13 0.5988176 0.5786799 1.0000000
#5001에서 10000명 데이터로 상관행렬 만들기
fdata2 <- data2 %>% mutate_all(as.numeric)
tc2 <- tetrachoric(fdata2)
tcc2 <-as.data.frame(tc2$rho)
tcc2
## B6_2 B6_3 B6_4 B6_6 B6_8 B6_9
## B6_2 1.0000000 0.6431286 0.6559586 0.6047201 0.6355097 0.5706985
## B6_3 0.6431286 1.0000000 0.7014715 0.6536939 0.7435952 0.4977546
## B6_4 0.6559586 0.7014715 1.0000000 0.5892777 0.6628593 0.6321918
## B6_6 0.6047201 0.6536939 0.5892777 1.0000000 0.7517755 0.5041825
## B6_8 0.6355097 0.7435952 0.6628593 0.7517755 1.0000000 0.5882458
## B6_9 0.5706985 0.4977546 0.6321918 0.5041825 0.5882458 1.0000000
## B6_10 0.4726389 0.3718790 0.4302827 0.4218484 0.4257507 0.4250811
## B6_12 0.5984993 0.7241564 0.6691890 0.6278776 0.7306937 0.5809721
## B6_14 0.5935565 0.7128325 0.6382914 0.6233900 0.7440889 0.5631836
## B6_15 0.4843246 0.6558430 0.5801444 0.5748999 0.7136327 0.5013333
## reB6_1 0.6543833 0.6615835 0.6492331 0.6283321 0.6730853 0.5286476
## reB6_5 0.6103634 0.6066529 0.6701260 0.5880978 0.6724016 0.5340590
## reB6_7 0.6232197 0.6116755 0.6807476 0.6409563 0.7360390 0.5668172
## reB6_11 0.5376545 0.5676068 0.5999101 0.5293429 0.6702760 0.5579655
## reB6_13 0.5747737 0.4435433 0.4891474 0.4304014 0.4603223 0.4563948
## B6_10 B6_12 B6_14 B6_15 reB6_1 reB6_5
## B6_2 0.4726389 0.5984993 0.5935565 0.4843246 0.6543833 0.6103634
## B6_3 0.3718790 0.7241564 0.7128325 0.6558430 0.6615835 0.6066529
## B6_4 0.4302827 0.6691890 0.6382914 0.5801444 0.6492331 0.6701260
## B6_6 0.4218484 0.6278776 0.6233900 0.5748999 0.6283321 0.5880978
## B6_8 0.4257507 0.7306937 0.7440889 0.7136327 0.6730853 0.6724016
## B6_9 0.4250811 0.5809721 0.5631836 0.5013333 0.5286476 0.5340590
## B6_10 1.0000000 0.4961247 0.4505362 0.4196250 0.4239548 0.4101019
## B6_12 0.4961247 1.0000000 0.8027057 0.6843897 0.5935654 0.5921330
## B6_14 0.4505362 0.8027057 1.0000000 0.8335362 0.6206567 0.5729236
## B6_15 0.4196250 0.6843897 0.8335362 1.0000000 0.6291964 0.5906804
## reB6_1 0.4239548 0.5935654 0.6206567 0.6291964 1.0000000 0.7447345
## reB6_5 0.4101019 0.5921330 0.5729236 0.5906804 0.7447345 1.0000000
## reB6_7 0.4508129 0.6210493 0.6205686 0.6360969 0.7803766 0.8667087
## reB6_11 0.4745637 0.6615553 0.5988085 0.6134621 0.6949073 0.7537720
## reB6_13 0.4854536 0.5100734 0.5101644 0.4795876 0.5455815 0.5835916
## reB6_7 reB6_11 reB6_13
## B6_2 0.6232197 0.5376545 0.5747737
## B6_3 0.6116755 0.5676068 0.4435433
## B6_4 0.6807476 0.5999101 0.4891474
## B6_6 0.6409563 0.5293429 0.4304014
## B6_8 0.7360390 0.6702760 0.4603223
## B6_9 0.5668172 0.5579655 0.4563948
## B6_10 0.4508129 0.4745637 0.4854536
## B6_12 0.6210493 0.6615553 0.5100734
## B6_14 0.6205686 0.5988085 0.5101644
## B6_15 0.6360969 0.6134621 0.4795876
## reB6_1 0.7803766 0.6949073 0.5455815
## reB6_5 0.8667087 0.7537720 0.5835916
## reB6_7 1.0000000 0.8170735 0.6098424
## reB6_11 0.8170735 1.0000000 0.5828405
## reB6_13 0.6098424 0.5828405 1.0000000
#plot
fa.parallel(fdata1)
## Parallel analysis suggests that the number of factors = 5 and the number of components = 2
fa.parallel(tcc1)
## Warning in fa.parallel(tcc1): It seems as if you are using a correlation
## matrix, but have not specified the number of cases. The number of subjects
## is arbitrarily set to be 100
## Parallel analysis suggests that the number of factors = 1 and the number of components = 1
scree(fdata1)
scree(tcc1)
efa <- fa(tcc1, 1, n.obs = 5000, fm="wls")
#EFA, 상관행렬을 넣어서 요인1개로 추정방법은 WLS = weighted least squares로 실행
efa
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = 1, n.obs = 5000, fm = "wls")
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS1 h2 u2 com
## B6_2 0.74 0.55 0.45 1
## B6_3 0.78 0.60 0.40 1
## B6_4 0.80 0.64 0.36 1
## B6_6 0.74 0.55 0.45 1
## B6_8 0.86 0.74 0.26 1
## B6_9 0.69 0.48 0.52 1
## B6_10 0.56 0.32 0.68 1
## B6_12 0.80 0.65 0.35 1
## B6_14 0.81 0.66 0.34 1
## B6_15 0.81 0.65 0.35 1
## reB6_1 0.77 0.60 0.40 1
## reB6_5 0.79 0.63 0.37 1
## reB6_7 0.82 0.68 0.32 1
## reB6_11 0.78 0.60 0.40 1
## reB6_13 0.66 0.43 0.57 1
##
## WLS1
## SS loadings 8.76
## Proportion Var 0.58
##
## Mean item complexity = 1
## Test of the hypothesis that 1 factor is sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 90 and the objective function was 2.74
##
## The root mean square of the residuals (RMSR) is 0.07
## The df corrected root mean square of the residuals is 0.07
##
## The harmonic number of observations is 5000 with the empirical chi square 4595.38 with prob < 0
## The total number of observations was 5000 with Likelihood Chi Square = 13656.7 with prob < 0
##
## Tucker Lewis Index of factoring reliability = 0.758
## RMSEA index = 0.174 and the 90 % confidence intervals are 0.171 0.176
## BIC = 12890.15
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy
## WLS1
## Correlation of (regression) scores with factors 0.98
## Multiple R square of scores with factors 0.96
## Minimum correlation of possible factor scores 0.92
print(efa$loadings, cutoff = 0.3)
##
## Loadings:
## WLS1
## B6_2 0.740
## B6_3 0.777
## B6_4 0.799
## B6_6 0.740
## B6_8 0.859
## B6_9 0.691
## B6_10 0.563
## B6_12 0.804
## B6_14 0.811
## B6_15 0.806
## reB6_1 0.772
## reB6_5 0.791
## reB6_7 0.823
## reB6_11 0.777
## reB6_13 0.658
##
## WLS1
## SS loadings 8.760
## Proportion Var 0.584
for(i in 1:6){
f <- fa(tcc1, i, n.obs = 5000, fm = "wls")
assign("f",i,f)
print(i)
print(f)
print(f$loadings, cutoff = 0.3)
}
## [1] 1
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = i, n.obs = 5000, fm = "wls")
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS1 h2 u2 com
## B6_2 0.74 0.55 0.45 1
## B6_3 0.78 0.60 0.40 1
## B6_4 0.80 0.64 0.36 1
## B6_6 0.74 0.55 0.45 1
## B6_8 0.86 0.74 0.26 1
## B6_9 0.69 0.48 0.52 1
## B6_10 0.56 0.32 0.68 1
## B6_12 0.80 0.65 0.35 1
## B6_14 0.81 0.66 0.34 1
## B6_15 0.81 0.65 0.35 1
## reB6_1 0.77 0.60 0.40 1
## reB6_5 0.79 0.63 0.37 1
## reB6_7 0.82 0.68 0.32 1
## reB6_11 0.78 0.60 0.40 1
## reB6_13 0.66 0.43 0.57 1
##
## WLS1
## SS loadings 8.76
## Proportion Var 0.58
##
## Mean item complexity = 1
## Test of the hypothesis that 1 factor is sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 90 and the objective function was 2.74
##
## The root mean square of the residuals (RMSR) is 0.07
## The df corrected root mean square of the residuals is 0.07
##
## The harmonic number of observations is 5000 with the empirical chi square 4595.38 with prob < 0
## The total number of observations was 5000 with Likelihood Chi Square = 13656.7 with prob < 0
##
## Tucker Lewis Index of factoring reliability = 0.758
## RMSEA index = 0.174 and the 90 % confidence intervals are 0.171 0.176
## BIC = 12890.15
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy
## WLS1
## Correlation of (regression) scores with factors 0.98
## Multiple R square of scores with factors 0.96
## Minimum correlation of possible factor scores 0.92
##
## Loadings:
## WLS1
## B6_2 0.740
## B6_3 0.777
## B6_4 0.799
## B6_6 0.740
## B6_8 0.859
## B6_9 0.691
## B6_10 0.563
## B6_12 0.804
## B6_14 0.811
## B6_15 0.806
## reB6_1 0.772
## reB6_5 0.791
## reB6_7 0.823
## reB6_11 0.777
## reB6_13 0.658
##
## WLS1
## SS loadings 8.760
## Proportion Var 0.584
## Loading required namespace: GPArotation
## [1] 2
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = i, n.obs = 5000, fm = "wls")
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS1 WLS2 h2 u2 com
## B6_2 0.51 0.27 0.54 0.458 1.5
## B6_3 0.81 0.00 0.66 0.341 1.0
## B6_4 0.58 0.26 0.63 0.365 1.4
## B6_6 0.64 0.14 0.56 0.442 1.1
## B6_8 0.73 0.17 0.75 0.250 1.1
## B6_9 0.48 0.25 0.47 0.526 1.5
## B6_10 0.51 0.07 0.32 0.675 1.0
## B6_12 0.87 -0.03 0.72 0.282 1.0
## B6_14 0.99 -0.14 0.79 0.213 1.0
## B6_15 0.80 0.05 0.69 0.306 1.0
## reB6_1 0.22 0.62 0.64 0.358 1.3
## reB6_5 -0.02 0.92 0.81 0.185 1.0
## reB6_7 -0.04 0.98 0.91 0.092 1.0
## reB6_11 0.16 0.70 0.68 0.319 1.1
## reB6_13 0.23 0.48 0.45 0.546 1.4
##
## WLS1 WLS2
## SS loadings 5.81 3.83
## Proportion Var 0.39 0.26
## Cumulative Var 0.39 0.64
## Proportion Explained 0.60 0.40
## Cumulative Proportion 0.60 1.00
##
## With factor correlations of
## WLS1 WLS2
## WLS1 1.00 0.74
## WLS2 0.74 1.00
##
## Mean item complexity = 1.2
## Test of the hypothesis that 2 factors are sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 76 and the objective function was 1.33
##
## The root mean square of the residuals (RMSR) is 0.04
## The df corrected root mean square of the residuals is 0.05
##
## The harmonic number of observations is 5000 with the empirical chi square 1777.08 with prob < 1.200085e-320
## The total number of observations was 5000 with Likelihood Chi Square = 6642.46 with prob < 0
##
## Tucker Lewis Index of factoring reliability = 0.861
## RMSEA index = 0.132 and the 90 % confidence intervals are 0.129 0.134
## BIC = 5995.16
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## WLS1 WLS2
## Correlation of (regression) scores with factors 0.97 0.98
## Multiple R square of scores with factors 0.95 0.96
## Minimum correlation of possible factor scores 0.90 0.91
##
## Loadings:
## WLS1 WLS2
## B6_2 0.510
## B6_3 0.814
## B6_4 0.585
## B6_6 0.639
## B6_8 0.733
## B6_9 0.479
## B6_10 0.514
## B6_12 0.866
## B6_14 0.988
## B6_15 0.796
## reB6_1 0.621
## reB6_5 0.917
## reB6_7 0.983
## reB6_11 0.702
## reB6_13 0.484
##
## WLS1 WLS2
## SS loadings 5.194 3.205
## Proportion Var 0.346 0.214
## Cumulative Var 0.346 0.560
## [1] 3
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = i, n.obs = 5000, fm = "wls")
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS1 WLS2 WLS3 h2 u2 com
## B6_2 0.21 0.19 0.50 0.62 0.384 1.7
## B6_3 0.73 0.02 0.09 0.66 0.338 1.0
## B6_4 0.39 0.23 0.30 0.64 0.360 2.5
## B6_6 0.56 0.15 0.11 0.56 0.445 1.2
## B6_8 0.72 0.20 0.01 0.77 0.227 1.2
## B6_9 0.28 0.21 0.32 0.49 0.508 2.7
## B6_10 0.18 -0.08 0.63 0.50 0.505 1.2
## B6_12 0.72 -0.02 0.20 0.72 0.283 1.2
## B6_14 0.99 -0.09 -0.03 0.83 0.174 1.0
## B6_15 0.84 0.11 -0.09 0.75 0.255 1.1
## reB6_1 0.18 0.61 0.09 0.65 0.353 1.2
## reB6_5 -0.03 0.91 0.05 0.83 0.169 1.0
## reB6_7 0.04 0.97 -0.07 0.93 0.074 1.0
## reB6_11 0.11 0.68 0.10 0.69 0.315 1.1
## reB6_13 -0.08 0.40 0.51 0.56 0.443 2.0
##
## WLS1 WLS2 WLS3
## SS loadings 4.69 3.66 1.82
## Proportion Var 0.31 0.24 0.12
## Cumulative Var 0.31 0.56 0.68
## Proportion Explained 0.46 0.36 0.18
## Cumulative Proportion 0.46 0.82 1.00
##
## With factor correlations of
## WLS1 WLS2 WLS3
## WLS1 1.00 0.70 0.63
## WLS2 0.70 1.00 0.56
## WLS3 0.63 0.56 1.00
##
## Mean item complexity = 1.4
## Test of the hypothesis that 3 factors are sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 63 and the objective function was 0.94
##
## The root mean square of the residuals (RMSR) is 0.03
## The df corrected root mean square of the residuals is 0.04
##
## The harmonic number of observations is 5000 with the empirical chi square 986.65 with prob < 5.6e-166
## The total number of observations was 5000 with Likelihood Chi Square = 4703.27 with prob < 0
##
## Tucker Lewis Index of factoring reliability = 0.882
## RMSEA index = 0.121 and the 90 % confidence intervals are 0.118 0.124
## BIC = 4166.69
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## WLS1 WLS2 WLS3
## Correlation of (regression) scores with factors 0.97 0.98 0.88
## Multiple R square of scores with factors 0.95 0.96 0.78
## Minimum correlation of possible factor scores 0.89 0.92 0.57
##
## Loadings:
## WLS1 WLS2 WLS3
## B6_2 0.496
## B6_3 0.735
## B6_4 0.391
## B6_6 0.556
## B6_8 0.719
## B6_9 0.319
## B6_10 0.625
## B6_12 0.721
## B6_14 0.987
## B6_15 0.840
## reB6_1 0.610
## reB6_5 0.907
## reB6_7 0.973
## reB6_11 0.684
## reB6_13 0.402 0.508
##
## WLS1 WLS2 WLS3
## SS loadings 3.926 2.994 1.182
## Proportion Var 0.262 0.200 0.079
## Cumulative Var 0.262 0.461 0.540
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.
## [1] 4
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = i, n.obs = 5000, fm = "wls")
##
## Warning: A Heywood case was detected.
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS1 WLS2 WLS3 WLS4 h2 u2 com
## B6_2 -0.02 0.07 0.79 0.09 0.71 0.289 1.0
## B6_3 0.51 0.01 0.40 -0.14 0.69 0.314 2.1
## B6_4 0.17 0.16 0.59 -0.04 0.69 0.306 1.3
## B6_6 0.37 0.13 0.37 -0.14 0.59 0.413 2.5
## B6_8 0.54 0.22 0.25 -0.15 0.78 0.215 2.0
## B6_9 0.14 0.16 0.47 0.04 0.51 0.495 1.4
## B6_10 0.24 -0.07 0.39 0.34 0.47 0.530 2.7
## B6_12 0.69 0.02 0.16 0.14 0.74 0.261 1.2
## B6_14 1.00 -0.03 -0.05 0.04 0.90 0.098 1.0
## B6_15 0.77 0.17 -0.02 -0.03 0.75 0.248 1.1
## reB6_1 0.06 0.56 0.27 -0.05 0.65 0.351 1.5
## reB6_5 -0.03 0.84 0.11 0.02 0.80 0.198 1.0
## reB6_7 0.00 1.02 -0.04 -0.05 0.98 0.017 1.0
## reB6_11 0.20 0.71 -0.06 0.18 0.73 0.268 1.3
## reB6_13 0.04 0.38 0.24 0.40 0.61 0.394 2.7
##
## WLS1 WLS2 WLS3 WLS4
## SS loadings 3.76 3.56 2.82 0.46
## Proportion Var 0.25 0.24 0.19 0.03
## Cumulative Var 0.25 0.49 0.68 0.71
## Proportion Explained 0.35 0.34 0.27 0.04
## Cumulative Proportion 0.35 0.69 0.96 1.00
##
## With factor correlations of
## WLS1 WLS2 WLS3 WLS4
## WLS1 1.00 0.65 0.69 0.14
## WLS2 0.65 1.00 0.66 0.13
## WLS3 0.69 0.66 1.00 0.17
## WLS4 0.14 0.13 0.17 1.00
##
## Mean item complexity = 1.6
## Test of the hypothesis that 4 factors are sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 51 and the objective function was 0.57
##
## The root mean square of the residuals (RMSR) is 0.02
## The df corrected root mean square of the residuals is 0.03
##
## The harmonic number of observations is 5000 with the empirical chi square 429.33 with prob < 2.9e-61
## The total number of observations was 5000 with Likelihood Chi Square = 2822.32 with prob < 0
##
## Tucker Lewis Index of factoring reliability = 0.913
## RMSEA index = 0.104 and the 90 % confidence intervals are 0.101 0.108
## BIC = 2387.95
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## WLS1 WLS2 WLS3 WLS4
## Correlation of (regression) scores with factors 0.97 1.00 0.94 0.74
## Multiple R square of scores with factors 0.95 0.99 0.88 0.54
## Minimum correlation of possible factor scores 0.90 0.99 0.76 0.09
##
## Loadings:
## WLS1 WLS2 WLS3 WLS4
## B6_2 0.790
## B6_3 0.510 0.398
## B6_4 0.594
## B6_6 0.365 0.374
## B6_8 0.538
## B6_9 0.472
## B6_10 0.394 0.339
## B6_12 0.693
## B6_14 1.002
## B6_15 0.768
## reB6_1 0.559
## reB6_5 0.842
## reB6_7 1.022
## reB6_11 0.712
## reB6_13 0.382 0.404
##
## WLS1 WLS2 WLS3 WLS4
## SS loadings 2.912 2.873 1.894 0.412
## Proportion Var 0.194 0.192 0.126 0.027
## Cumulative Var 0.194 0.386 0.512 0.539
## [1] 5
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = i, n.obs = 5000, fm = "wls")
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS2 WLS1 WLS3 WLS4 WLS5 h2 u2 com
## B6_2 0.10 -0.05 0.34 0.48 0.18 0.69 0.309 2.3
## B6_3 0.03 0.50 -0.08 0.42 0.10 0.73 0.275 2.1
## B6_4 0.17 0.17 0.11 0.49 0.12 0.71 0.290 1.8
## B6_6 0.12 0.12 0.11 0.12 0.53 0.65 0.347 1.4
## B6_8 0.20 0.30 0.04 0.03 0.55 0.88 0.124 1.9
## B6_9 0.17 0.10 0.21 0.30 0.15 0.50 0.496 3.3
## B6_10 -0.10 0.12 0.60 0.04 0.16 0.52 0.482 1.3
## B6_12 0.02 0.64 0.20 0.11 0.03 0.73 0.267 1.3
## B6_14 -0.03 0.94 0.03 -0.02 0.04 0.91 0.094 1.0
## B6_15 0.16 0.72 -0.03 0.00 0.08 0.76 0.239 1.1
## reB6_1 0.56 0.11 -0.02 0.30 -0.01 0.67 0.326 1.6
## reB6_5 0.82 0.01 0.03 0.12 -0.02 0.81 0.191 1.0
## reB6_7 0.97 -0.04 -0.01 -0.05 0.15 0.99 0.014 1.1
## reB6_11 0.67 0.25 0.14 -0.02 -0.12 0.72 0.275 1.4
## reB6_13 0.36 0.07 0.53 0.03 -0.10 0.62 0.380 1.9
##
## WLS2 WLS1 WLS3 WLS4 WLS5
## SS loadings 3.41 3.18 1.41 1.58 1.31
## Proportion Var 0.23 0.21 0.09 0.11 0.09
## Cumulative Var 0.23 0.44 0.53 0.64 0.73
## Proportion Explained 0.31 0.29 0.13 0.15 0.12
## Cumulative Proportion 0.31 0.60 0.73 0.88 1.00
##
## With factor correlations of
## WLS2 WLS1 WLS3 WLS4 WLS5
## WLS2 1.00 0.60 0.49 0.49 0.50
## WLS1 0.60 1.00 0.51 0.53 0.62
## WLS3 0.49 0.51 1.00 0.43 0.31
## WLS4 0.49 0.53 0.43 1.00 0.44
## WLS5 0.50 0.62 0.31 0.44 1.00
##
## Mean item complexity = 1.6
## Test of the hypothesis that 5 factors are sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 40 and the objective function was 0.33
##
## The root mean square of the residuals (RMSR) is 0.01
## The df corrected root mean square of the residuals is 0.02
##
## The harmonic number of observations is 5000 with the empirical chi square 232.56 with prob < 5.4e-29
## The total number of observations was 5000 with Likelihood Chi Square = 1640.43 with prob < 1.200001e-318
##
## Tucker Lewis Index of factoring reliability = 0.936
## RMSEA index = 0.09 and the 90 % confidence intervals are 0.086 0.093
## BIC = 1299.74
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## WLS2 WLS1 WLS3 WLS4 WLS5
## Correlation of (regression) scores with factors 0.99 0.97 0.85 0.88 0.90
## Multiple R square of scores with factors 0.98 0.94 0.73 0.77 0.81
## Minimum correlation of possible factor scores 0.97 0.89 0.45 0.54 0.62
##
## Loadings:
## WLS2 WLS1 WLS3 WLS4 WLS5
## B6_2 0.341 0.481
## B6_3 0.500 0.425
## B6_4 0.487
## B6_6 0.526
## B6_8 0.547
## B6_9
## B6_10 0.605
## B6_12 0.639
## B6_14 0.940
## B6_15 0.725
## reB6_1 0.563 0.305
## reB6_5 0.817
## reB6_7 0.965
## reB6_11 0.667
## reB6_13 0.358 0.527
##
## WLS2 WLS1 WLS3 WLS4 WLS5
## SS loadings 2.649 2.305 0.899 0.880 0.740
## Proportion Var 0.177 0.154 0.060 0.059 0.049
## Cumulative Var 0.177 0.330 0.390 0.449 0.498
## Warning in GPFoblq(L, Tmat = Tmat, normalize = normalize, eps = eps, maxit
## = maxit, : convergence not obtained in GPFoblq. 1000 iterations used.
## [1] 6
## Factor Analysis using method = wls
## Call: fa(r = tcc1, nfactors = i, n.obs = 5000, fm = "wls")
## Standardized loadings (pattern matrix) based upon correlation matrix
## WLS2 WLS1 WLS5 WLS3 WLS4 WLS6 h2 u2 com
## B6_2 0.04 0.05 0.15 0.43 0.41 -0.17 0.73 0.271 2.6
## B6_3 0.03 0.38 0.13 0.47 -0.07 0.02 0.71 0.289 2.2
## B6_4 0.17 0.00 0.14 0.61 0.07 0.09 0.75 0.252 1.3
## B6_6 0.06 0.10 0.57 0.12 0.10 -0.05 0.64 0.362 1.3
## B6_8 0.13 0.14 0.73 0.04 0.01 0.04 0.91 0.090 1.2
## B6_9 0.15 0.01 0.19 0.33 0.20 0.05 0.51 0.495 2.9
## B6_10 -0.11 0.05 0.24 0.05 0.54 0.12 0.50 0.505 1.6
## B6_12 0.03 0.39 0.12 0.25 0.12 0.32 0.80 0.205 3.2
## B6_14 -0.03 0.97 0.02 -0.02 0.03 0.01 0.96 0.044 1.0
## B6_15 0.16 0.69 0.10 0.02 -0.01 -0.02 0.76 0.244 1.2
## reB6_1 0.54 0.17 -0.02 0.27 0.02 -0.11 0.67 0.327 1.8
## reB6_5 0.81 0.06 -0.03 0.10 0.05 -0.06 0.82 0.177 1.1
## reB6_7 0.90 -0.03 0.20 -0.05 -0.01 -0.02 0.96 0.039 1.1
## reB6_11 0.70 0.08 -0.04 0.05 0.09 0.25 0.77 0.226 1.3
## reB6_13 0.34 0.13 -0.09 -0.03 0.56 0.02 0.64 0.362 1.8
##
## WLS2 WLS1 WLS5 WLS3 WLS4 WLS6
## SS loadings 3.20 2.58 1.80 1.85 1.35 0.33
## Proportion Var 0.21 0.17 0.12 0.12 0.09 0.02
## Cumulative Var 0.21 0.39 0.51 0.63 0.72 0.74
## Proportion Explained 0.29 0.23 0.16 0.17 0.12 0.03
## Cumulative Proportion 0.29 0.52 0.68 0.85 0.97 1.00
##
## With factor correlations of
## WLS2 WLS1 WLS5 WLS3 WLS4 WLS6
## WLS2 1.00 0.56 0.59 0.52 0.49 0.13
## WLS1 0.56 1.00 0.72 0.60 0.46 0.34
## WLS5 0.59 0.72 1.00 0.58 0.38 0.15
## WLS3 0.52 0.60 0.58 1.00 0.48 0.08
## WLS4 0.49 0.46 0.38 0.48 1.00 0.19
## WLS6 0.13 0.34 0.15 0.08 0.19 1.00
##
## Mean item complexity = 1.7
## Test of the hypothesis that 6 factors are sufficient.
##
## The degrees of freedom for the null model are 105 and the objective function was 13.14 with Chi Square of 65621.3
## The degrees of freedom for the model are 30 and the objective function was 0.21
##
## The root mean square of the residuals (RMSR) is 0.01
## The df corrected root mean square of the residuals is 0.02
##
## The harmonic number of observations is 5000 with the empirical chi square 132.66 with prob < 7.2e-15
## The total number of observations was 5000 with Likelihood Chi Square = 1040.01 with prob < 1.8e-199
##
## Tucker Lewis Index of factoring reliability = 0.946
## RMSEA index = 0.082 and the 90 % confidence intervals are 0.078 0.086
## BIC = 784.49
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## WLS2 WLS1 WLS5 WLS3 WLS4
## Correlation of (regression) scores with factors 0.98 0.98 0.95 0.90 0.85
## Multiple R square of scores with factors 0.97 0.96 0.90 0.82 0.73
## Minimum correlation of possible factor scores 0.93 0.93 0.79 0.63 0.46
## WLS6
## Correlation of (regression) scores with factors 0.73
## Multiple R square of scores with factors 0.54
## Minimum correlation of possible factor scores 0.07
##
## Loadings:
## WLS2 WLS1 WLS5 WLS3 WLS4 WLS6
## B6_2 0.426 0.407
## B6_3 0.379 0.466
## B6_4 0.609
## B6_6 0.567
## B6_8 0.727
## B6_9 0.328
## B6_10 0.543
## B6_12 0.386 0.324
## B6_14 0.975
## B6_15 0.691
## reB6_1 0.541
## reB6_5 0.808
## reB6_7 0.900
## reB6_11 0.705
## reB6_13 0.337 0.558
##
## WLS2 WLS1 WLS5 WLS3 WLS4 WLS6
## SS loadings 2.480 1.808 1.076 1.043 0.858 0.244
## Proportion Var 0.165 0.121 0.072 0.070 0.057 0.016
## Cumulative Var 0.165 0.286 0.358 0.427 0.484 0.501
4. CFA
library(lavaan)
## Warning: package 'lavaan' was built under R version 3.5.3
## This is lavaan 0.6-5
## lavaan is BETA software! Please report any bugs.
##
## Attaching package: 'lavaan'
## The following object is masked from 'package:psych':
##
## cor2cov
sd <- as.data.frame(apply(data1, 2, sd))
nsd <- sd$`apply(data1, 2, sd)`
covdata <- cor2cov(tc1$rho, nsd)
model1 <- '
Depression =~ NA*B6_2 + B6_3 + B6_4 + B6_6 + B6_8 + B6_9 + B6_10 + B6_12 + B6_14 + B6_15 + reB6_1 + reB6_5 + reB6_7 + reB6_11 + reB6_13
'
#model syntax = x1 ~ m1*1 # Means / x1 ~~ r1*x1 #Variance / x1 ~~ 0*x2 #residual covariate
#NA*는 첫번째 indicator를 free로 추정하는 방법.
# std.lv = TRUE 는 요인 분산 1로 고정하는 명령어.
fit1 <- cfa(model1, sample.cov = covdata, sample.nobs = 5000, std.lv = T, mimic = "Mplus")
summary(fit1, fit.measures = TRUE)
## lavaan 0.6-5 ended normally after 64 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free parameters 45
##
## Number of observations 5000
##
## Model Test User Model:
##
## Test statistic 13663.220
## Degrees of freedom 90
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 65711.103
## Degrees of freedom 105
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.793
## Tucker-Lewis Index (TLI) 0.759
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15221.011
## Loglikelihood unrestricted model (H1) -8389.401
##
## Akaike (AIC) 30532.023
## Bayesian (BIC) 30825.297
## Sample-size adjusted Bayesian (BIC) 30682.302
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.174
## 90 Percent confidence interval - lower 0.171
## 90 Percent confidence interval - upper 0.176
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.059
##
## Parameter Estimates:
##
## Information Observed
## Observed information based on Hessian
## Standard errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## Depression =~
## B6_2 0.358 0.006 59.204 0.000
## B6_3 0.279 0.004 65.412 0.000
## B6_4 0.370 0.006 66.856 0.000
## B6_6 0.269 0.004 61.158 0.000
## B6_8 0.291 0.004 76.154 0.000
## B6_9 0.242 0.004 54.666 0.000
## B6_10 0.241 0.006 41.986 0.000
## B6_12 0.307 0.005 68.133 0.000
## B6_14 0.317 0.005 69.608 0.000
## B6_15 0.296 0.004 69.361 0.000
## reB6_1 0.371 0.006 63.980 0.000
## reB6_5 0.388 0.006 66.046 0.000
## reB6_7 0.397 0.006 70.398 0.000
## reB6_11 0.368 0.006 64.589 0.000
## reB6_13 0.323 0.006 50.672 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .B6_2 0.000 0.007 0.000 1.000
## .B6_3 0.000 0.005 0.000 1.000
## .B6_4 0.000 0.007 0.000 1.000
## .B6_6 0.000 0.005 0.000 1.000
## .B6_8 0.000 0.005 0.000 1.000
## .B6_9 0.000 0.005 0.000 1.000
## .B6_10 0.000 0.006 0.000 1.000
## .B6_12 0.000 0.005 0.000 1.000
## .B6_14 0.000 0.005 0.000 1.000
## .B6_15 0.000 0.005 0.000 1.000
## .reB6_1 0.000 0.007 0.000 1.000
## .reB6_5 0.000 0.007 0.000 1.000
## .reB6_7 0.000 0.007 0.000 1.000
## .reB6_11 0.000 0.007 0.000 1.000
## .reB6_13 0.000 0.007 0.000 1.000
## Depression 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .B6_2 0.113 0.002 47.392 0.000
## .B6_3 0.049 0.001 46.149 0.000
## .B6_4 0.081 0.002 46.162 0.000
## .B6_6 0.058 0.001 47.129 0.000
## .B6_8 0.029 0.001 43.197 0.000
## .B6_9 0.065 0.001 48.029 0.000
## .B6_10 0.129 0.003 48.990 0.000
## .B6_12 0.051 0.001 45.345 0.000
## .B6_14 0.050 0.001 44.280 0.000
## .B6_15 0.045 0.001 45.072 0.000
## .reB6_1 0.094 0.002 46.346 0.000
## .reB6_5 0.091 0.002 45.043 0.000
## .reB6_7 0.075 0.002 43.707 0.000
## .reB6_11 0.089 0.002 46.009 0.000
## .reB6_13 0.144 0.003 48.278 0.000
## Depression 1.000
standardizedSolution(fit1)
## lhs op rhs est.std se z pvalue ci.lower ci.upper
## 1 Depression =~ B6_2 0.729 0.007 104.518 0 0.715 0.743
## 2 Depression =~ B6_3 0.782 0.006 133.465 0 0.771 0.794
## 3 Depression =~ B6_4 0.793 0.006 141.959 0 0.782 0.804
## 4 Depression =~ B6_6 0.746 0.007 112.750 0 0.733 0.759
## 5 Depression =~ B6_8 0.863 0.004 215.995 0 0.855 0.871
## 6 Depression =~ B6_9 0.687 0.008 88.064 0 0.671 0.702
## 7 Depression =~ B6_10 0.556 0.010 55.121 0 0.536 0.576
## 8 Depression =~ B6_12 0.805 0.005 149.293 0 0.794 0.815
## 9 Depression =~ B6_14 0.817 0.005 158.293 0 0.807 0.827
## 10 Depression =~ B6_15 0.814 0.005 157.402 0 0.804 0.824
## 11 Depression =~ reB6_1 0.771 0.006 125.875 0 0.759 0.783
## 12 Depression =~ reB6_5 0.790 0.006 136.151 0 0.778 0.801
## 13 Depression =~ reB6_7 0.824 0.005 163.415 0 0.814 0.833
## 14 Depression =~ reB6_11 0.776 0.006 128.835 0 0.765 0.788
## 15 Depression =~ reB6_13 0.648 0.009 75.906 0 0.632 0.665
## 16 B6_2 ~~ B6_2 0.469 0.010 46.116 0 0.449 0.489
## 17 B6_3 ~~ B6_3 0.388 0.009 42.287 0 0.370 0.406
## 18 B6_4 ~~ B6_4 0.371 0.009 41.781 0 0.353 0.388
## 19 B6_6 ~~ B6_6 0.443 0.010 44.927 0 0.424 0.463
## 20 B6_8 ~~ B6_8 0.255 0.007 36.963 0 0.241 0.269
## 21 B6_9 ~~ B6_9 0.528 0.011 49.340 0 0.507 0.549
## 22 B6_10 ~~ B6_10 0.691 0.011 61.528 0 0.669 0.713
## 23 B6_12 ~~ B6_12 0.353 0.009 40.685 0 0.336 0.370
## 24 B6_14 ~~ B6_14 0.332 0.008 39.402 0 0.316 0.349
## 25 B6_15 ~~ B6_15 0.337 0.008 40.087 0 0.321 0.354
## 26 reB6_1 ~~ reB6_1 0.406 0.009 42.999 0 0.387 0.424
## 27 reB6_5 ~~ reB6_5 0.377 0.009 41.111 0 0.359 0.394
## 28 reB6_7 ~~ reB6_7 0.322 0.008 38.745 0 0.305 0.338
## 29 reB6_11 ~~ reB6_11 0.397 0.009 42.464 0 0.379 0.416
## 30 reB6_13 ~~ reB6_13 0.580 0.011 52.313 0 0.558 0.601
## 31 Depression ~~ Depression 1.000 0.000 NA NA 1.000 1.000
## 32 B6_2 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 33 B6_3 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 34 B6_4 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 35 B6_6 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 36 B6_8 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 37 B6_9 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 38 B6_10 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 39 B6_12 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 40 B6_14 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 41 B6_15 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 42 reB6_1 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 43 reB6_5 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 44 reB6_7 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 45 reB6_11 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 46 reB6_13 ~1 0.000 0.014 0.000 1 -0.028 0.028
## 47 Depression ~1 0.000 0.000 NA NA 0.000 0.000
parameterEstimates(fit1,standardized = TRUE)
## lhs op rhs est se z pvalue ci.lower ci.upper
## 1 Depression =~ B6_2 0.358 0.006 59.204 0 0.347 0.370
## 2 Depression =~ B6_3 0.279 0.004 65.412 0 0.270 0.287
## 3 Depression =~ B6_4 0.370 0.006 66.856 0 0.359 0.381
## 4 Depression =~ B6_6 0.269 0.004 61.158 0 0.260 0.277
## 5 Depression =~ B6_8 0.291 0.004 76.154 0 0.284 0.299
## 6 Depression =~ B6_9 0.242 0.004 54.666 0 0.233 0.250
## 7 Depression =~ B6_10 0.241 0.006 41.986 0 0.229 0.252
## 8 Depression =~ B6_12 0.307 0.005 68.133 0 0.298 0.316
## 9 Depression =~ B6_14 0.317 0.005 69.608 0 0.308 0.326
## 10 Depression =~ B6_15 0.296 0.004 69.361 0 0.287 0.304
## 11 Depression =~ reB6_1 0.371 0.006 63.980 0 0.360 0.382
## 12 Depression =~ reB6_5 0.388 0.006 66.046 0 0.376 0.399
## 13 Depression =~ reB6_7 0.397 0.006 70.398 0 0.386 0.408
## 14 Depression =~ reB6_11 0.368 0.006 64.589 0 0.357 0.379
## 15 Depression =~ reB6_13 0.323 0.006 50.672 0 0.310 0.335
## 16 B6_2 ~~ B6_2 0.113 0.002 47.392 0 0.109 0.118
## 17 B6_3 ~~ B6_3 0.049 0.001 46.149 0 0.047 0.051
## 18 B6_4 ~~ B6_4 0.081 0.002 46.162 0 0.077 0.084
## 19 B6_6 ~~ B6_6 0.058 0.001 47.129 0 0.055 0.060
## 20 B6_8 ~~ B6_8 0.029 0.001 43.197 0 0.028 0.030
## 21 B6_9 ~~ B6_9 0.065 0.001 48.029 0 0.063 0.068
## 22 B6_10 ~~ B6_10 0.129 0.003 48.990 0 0.124 0.134
## 23 B6_12 ~~ B6_12 0.051 0.001 45.345 0 0.049 0.054
## 24 B6_14 ~~ B6_14 0.050 0.001 44.280 0 0.048 0.052
## 25 B6_15 ~~ B6_15 0.045 0.001 45.072 0 0.043 0.047
## 26 reB6_1 ~~ reB6_1 0.094 0.002 46.346 0 0.090 0.098
## 27 reB6_5 ~~ reB6_5 0.091 0.002 45.043 0 0.087 0.095
## 28 reB6_7 ~~ reB6_7 0.075 0.002 43.707 0 0.071 0.078
## 29 reB6_11 ~~ reB6_11 0.089 0.002 46.009 0 0.085 0.093
## 30 reB6_13 ~~ reB6_13 0.144 0.003 48.278 0 0.138 0.150
## 31 Depression ~~ Depression 1.000 0.000 NA NA 1.000 1.000
## 32 B6_2 ~1 0.000 0.007 0.000 1 -0.014 0.014
## 33 B6_3 ~1 0.000 0.005 0.000 1 -0.010 0.010
## 34 B6_4 ~1 0.000 0.007 0.000 1 -0.013 0.013
## 35 B6_6 ~1 0.000 0.005 0.000 1 -0.010 0.010
## 36 B6_8 ~1 0.000 0.005 0.000 1 -0.009 0.009
## 37 B6_9 ~1 0.000 0.005 0.000 1 -0.010 0.010
## 38 B6_10 ~1 0.000 0.006 0.000 1 -0.012 0.012
## 39 B6_12 ~1 0.000 0.005 0.000 1 -0.011 0.011
## 40 B6_14 ~1 0.000 0.005 0.000 1 -0.011 0.011
## 41 B6_15 ~1 0.000 0.005 0.000 1 -0.010 0.010
## 42 reB6_1 ~1 0.000 0.007 0.000 1 -0.013 0.013
## 43 reB6_5 ~1 0.000 0.007 0.000 1 -0.014 0.014
## 44 reB6_7 ~1 0.000 0.007 0.000 1 -0.013 0.013
## 45 reB6_11 ~1 0.000 0.007 0.000 1 -0.013 0.013
## 46 reB6_13 ~1 0.000 0.007 0.000 1 -0.014 0.014
## 47 Depression ~1 0.000 0.000 NA NA 0.000 0.000
## std.lv std.all std.nox
## 1 0.358 0.729 0.729
## 2 0.279 0.782 0.782
## 3 0.370 0.793 0.793
## 4 0.269 0.746 0.746
## 5 0.291 0.863 0.863
## 6 0.242 0.687 0.687
## 7 0.241 0.556 0.556
## 8 0.307 0.805 0.805
## 9 0.317 0.817 0.817
## 10 0.296 0.814 0.814
## 11 0.371 0.771 0.771
## 12 0.388 0.790 0.790
## 13 0.397 0.824 0.824
## 14 0.368 0.776 0.776
## 15 0.323 0.648 0.648
## 16 0.113 0.469 0.469
## 17 0.049 0.388 0.388
## 18 0.081 0.371 0.371
## 19 0.058 0.443 0.443
## 20 0.029 0.255 0.255
## 21 0.065 0.528 0.528
## 22 0.129 0.691 0.691
## 23 0.051 0.353 0.353
## 24 0.050 0.332 0.332
## 25 0.045 0.337 0.337
## 26 0.094 0.406 0.406
## 27 0.091 0.377 0.377
## 28 0.075 0.322 0.322
## 29 0.089 0.397 0.397
## 30 0.144 0.580 0.580
## 31 1.000 1.000 1.000
## 32 0.000 0.000 0.000
## 33 0.000 0.000 0.000
## 34 0.000 0.000 0.000
## 35 0.000 0.000 0.000
## 36 0.000 0.000 0.000
## 37 0.000 0.000 0.000
## 38 0.000 0.000 0.000
## 39 0.000 0.000 0.000
## 40 0.000 0.000 0.000
## 41 0.000 0.000 0.000
## 42 0.000 0.000 0.000
## 43 0.000 0.000 0.000
## 44 0.000 0.000 0.000
## 45 0.000 0.000 0.000
## 46 0.000 0.000 0.000
## 47 0.000 0.000 0.000
library(semPlot)
## Warning: package 'semPlot' was built under R version 3.5.3
semPaths(fit1,
what = "est", style = "lisrel", layout = "tree", intStyle = "multi",
intercepts = FALSE, residuals = TRUE, esize = 1, asize = 3,
edge.color = "black", edge.label.cex = 0.9, nDigits = 2,
sizeMan = 4, sizeMan2 = 5, nCharNodes = 0)
library(MplusAutomation)
## Version: 0.7-2
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 1-18. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").
efa1 <- readModels("C:/Users/LG/Desktop/efa/efa1.out")
## Reading model: C:/Users/LG/Desktop/efa/efa1.out
## Error extracting model summaries in output file: C:/Users/LG/Desktop/efa/efa1.out
## <simpleError in extractSummaries_1file(outfiletext, curfile, input = inp): Unable to locate section headers for EFA model fit statistics>
## Error extracting latent class counts in output file: C:/Users/LG/Desktop/efa/efa1.out
## <simpleError in if (missing(summaries) || summaries$NCategoricalLatentVars == 1 || is.na(summaries$NCategoricalLatentVars)) { if (missing(summaries) || is.null(summaries$Mplus.version) || as.numeric(summaries$Mplus.version) < 7.3) { modelCounts <- getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASSES$", outfiletext) ppCounts <- getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASS PATTERNS$", outfiletext) mostLikelyCounts <- getSection("^CLASSIFICATION OF INDIVIDUALS BASED ON THEIR MOST LIKELY LATENT CLASS MEMBERSHIP$", outfiletext) } else { modelCounts <- getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASSES$::^BASED ON THE ESTIMATED MODEL$", outfiletext) ppCounts <- getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASSES$::^BASED ON ESTIMATED POSTERIOR PROBABILITIES$", outfiletext) mostLikelyCounts <- getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASSES$::^BASED ON THEIR MOST LIKELY LATENT CLASS MEMBERSHIP$", outfiletext) } countlist[["modelEstimated"]] <- getClassCols(modelCounts) countlist[["posteriorProb"]] <- getClassCols(ppCounts) countlist[["mostLikely"]] <- getClassCols(mostLikelyCounts) mostLikelyProbs <- getSection("^Average Latent Class Probabilities for Most Likely Latent Class Membership \\((Row|Column)\\)$", outfiletext) if (length(mostLikelyProbs) > 1L) { mostLikelyProbs <- mostLikelyProbs[-1L] } countlist[["avgProbs.mostLikely"]] <- unlabeledMatrixExtract(mostLikelyProbs, filename) classificationProbs <- getSection("^Classification Probabilities for the Most Likely Latent Class Membership \\((Column|Row)\\)$", outfiletext) if (length(classificationProbs) > 1L) { classificationProbs <- classificationProbs[-1L] } countlist[["classificationProbs.mostLikely"]] <- unlabeledMatrixExtract(classificationProbs, filename) classificationLogitProbs <- getSection("^Logits for the Classification Probabilities for the Most Likely Latent Class Membership \\((Column|Row)\\)$", outfiletext) if (length(classificationLogitProbs) > 1L) { classificationLogitProbs <- classificationLogitProbs[-1L] } countlist[["logitProbs.mostLikely"]] <- unlabeledMatrixExtract(classificationLogitProbs, filename)} else { getClassCols_lta <- function(sectiontext) { numberLines <- grep("^\\s*([a-zA-Z0-9]+)?(\\s+[0-9\\.-]{1,}){1,}$", sectiontext, perl = TRUE) if (length(numberLines) > 0) { parsedlines <- strsplit(trimSpace(sectiontext[numberLines]), "\\s+", perl = TRUE) num_values <- sapply(parsedlines, length) if (length(unique(num_values)) == 1) { counts <- data.frame(t(sapply(parsedlines, as.numeric)), stringsAsFactors = FALSE) } else { parsedlines[which(num_values != max(num_values))] <- lapply(parsedlines[which(num_values != max(num_values))], function(x) { c(rep(NA, (max(num_values) - length(x))), x) }) counts <- do.call(rbind, parsedlines) counts[, 1] <- inverse.rle(list(lengths = diff(c(which(!is.na(counts[, 1])), (nrow(counts) + 1))), values = counts[, 1][complete.cases(counts[, 1])])) counts <- data.frame(counts, stringsAsFactors = FALSE) counts[, 2:4] <- lapply(counts[, 2:4], as.numeric) } return(counts) } else { return(NULL) } } if (is.null(summaries$Mplus.version) || as.numeric(summaries$Mplus.version) < 7.3) { warning("MplusAutomation currently only reads output with multiple categorical latent variables for Mplus version 8 and up.") } countlist[["modelEstimated"]] <- getClassCols_lta(getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR EACH LATENT CLASS VARIABLE$::^BASED ON THE ESTIMATED MODEL$", outfiletext)) countlist[["posteriorProb"]] <- getClassCols_lta(getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR EACH LATENT CLASS VARIABLE$::^BASED ON ESTIMATED POSTERIOR PROBABILITIES$", outfiletext)) countlist[["mostLikely"]] <- getClassCols_lta(getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR EACH LATENT CLASS VARIABLE$::^BASED ON THEIR MOST LIKELY LATENT CLASS PATTERN$", outfiletext)) countlist[c("modelEstimated", "posteriorProb", "mostLikely")] <- lapply(countlist[c("modelEstimated", "posteriorProb", "mostLikely")], setNames, c("variable", "class", "count", "proportion")) countlist[["modelEstimated.patterns"]] <- getClassCols_lta(getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASS PATTERNS$::^BASED ON THE ESTIMATED MODEL$", outfiletext)) countlist[["posteriorProb.patterns"]] <- getClassCols_lta(getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASS PATTERNS$::^BASED ON ESTIMATED POSTERIOR PROBABILITIES$", outfiletext)) countlist[["mostLikely.patterns"]] <- getClassCols_lta(getSection("^FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASS PATTERNS$::^BASED ON THEIR MOST LIKELY LATENT CLASS PATTERN$", outfiletext)) countlist[c("modelEstimated.patterns", "posteriorProb.patterns", "mostLikely.patterns")] <- lapply(countlist[c("modelEstimated.patterns", "posteriorProb.patterns", "mostLikely.patterns")], setNames, c(paste0("class.", unique(countlist[["modelEstimated"]]$variable)), "count", "proportion")) avgProbs <- outfiletext[(grep("^Average Latent Class Probabilities for Most Likely Latent Class Pattern \\((Row|Column)\\)$", outfiletext) + 2):grep("^MODEL RESULTS$", outfiletext) - 1] column_headers <- strsplit(trimws(grep("\\s*Latent Class\\s{2,}", avgProbs, value = TRUE)), "\\s+", perl = TRUE)[[1]][-1] variable_pattern_rows <- grep(paste(c("^(\\s{2,}\\d+){", length(column_headers), "}$"), collapse = ""), avgProbs, perl = TRUE) variable_pattern_rows <- variable_pattern_rows[!c(FALSE, diff(variable_pattern_rows) != 1)] variable_patterns <- avgProbs[variable_pattern_rows] variable_patterns <- data.frame(t(sapply(strsplit(trimws(variable_patterns), "\\s+", perl = TRUE), as.numeric))) names(variable_patterns) <- c("Latent Class Pattern No.", column_headers[-1]) probs <- grep(paste(c("^\\s+\\d{1,}(\\s{2,}[0-9\\.-]+)+$"), collapse = ""), avgProbs[(variable_pattern_rows[length(variable_pattern_rows)] + 1):length(avgProbs)], perl = TRUE, value = TRUE) if (length(probs)%%nrow(variable_patterns) > 1) { for (i in 2:(length(probs)%%nrow(variable_patterns))) { probs[1:(nrow(variable_patterns) + 1)] <- paste(probs[1:(nrow(variable_patterns) + 1)], substring(probs[((i - 1) * (nrow(variable_patterns) + 1) + 1):(i * (nrow(variable_patterns) + 1))], first = 8)) } probs <- probs[1:nrow(variable_patterns)] } probs <- t(sapply(strsplit(trimws(probs[-1]), "\\s+", perl = TRUE), as.numeric))[, -1] countlist[["avgProbs.mostLikely"]] <- probs countlist[["avgProbs.mostLikely.patterns"]] <- variable_patterns countlist[["classificationProbs.mostLikely"]] <- NULL countlist[["logitProbs.mostLikely"]] <- NULL transitionProbs <- getSection("^LATENT TRANSITION PROBABILITIES BASED ON THE ESTIMATED MODEL$", outfiletext) section_starts <- grep("\\(Columns\\)$", transitionProbs) transitionProbs <- mapply(FUN = function(begin, end) { probs <- grep("^\\s+\\d{1,}(\\s{2,}[0-9\\.-]{2,}){1,}$", transitionProbs[begin:end], perl = TRUE, value = TRUE) probs <- do.call(rbind, strsplit(trimws(probs), "\\s+", perl = TRUE))[, -1] cbind(paste(gsub("\\s+(\\w+) Classes.*$", "\\1", transitionProbs[begin]), ".", rep(c(1:nrow(probs)), ncol(probs)), sep = ""), paste(gsub(".+?by (\\w+) Classes.*$", "\\1", transitionProbs[begin]), ".", as.vector(sapply(1:ncol(probs), rep, nrow(probs))), sep = ""), as.vector(probs)) }, begin = section_starts, end = c(section_starts[-1], length(transitionProbs)), SIMPLIFY = FALSE) if (length(transitionProbs) > 1) { transitionProbs <- do.call(rbind, transitionProbs) } else { transitionProbs <- transitionProbs[[1]] } transitionProbs <- data.frame(transitionProbs, stringsAsFactors = FALSE) names(transitionProbs) <- c("from", "to", "probability") transitionProbs$probability <- as.numeric(transitionProbs$probability) countlist[["transitionProbs"]] <- transitionProbs}: TRUE/FALSE가 필요한 곳에 값이 없습니다>
print(efa1$summaries)
## list()
5. Regression (moderation)
setwd("C:/Users/LG/Desktop/efa")
library(readxl)
rdata <- read_excel(path="regre.xlsx", sheet="regre", col_names=TRUE)
head(rdata)
## # A tibble: 6 x 49
## B6_1 B6_2 B6_3 B6_4 B6_5 B6_6 B6_7 B6_8 B6_9 B6_10 B6_11 B6_12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 1 2 1 2 2 2 1 2 1 2 1
## 2 1 1 2 2 1 2 1 2 2 2 1 1
## 3 1 1 2 2 1 2 1 2 2 2 1 1
## 4 2 2 1 2 2 1 2 2 2 2 2 2
## 5 1 1 2 1 1 2 1 2 2 2 1 2
## 6 1 1 2 2 1 2 1 2 2 1 1 2
## # ... with 37 more variables: B6_13 <dbl>, B6_14 <dbl>, B6_15 <dbl>,
## # C5_1 <dbl>, C5_2 <dbl>, C5_3_1 <dbl>, C5_3_2 <dbl>, C5_3_3 <dbl>,
## # C5_4 <dbl>, C5_5 <dbl>, C5_6 <dbl>, C5_7 <dbl>, C5_8 <dbl>,
## # C5_9 <dbl>, C5_10 <dbl>, 공변량 <lgl>, ID_A3 <dbl>, ID_A4_2 <dbl>,
## # ID_A5 <dbl>, ID_A6_1 <dbl>, ID_A7 <dbl>, B1 <dbl>, VAR00003 <lgl>,
## # wg_h <dbl>, ws_h <dbl>, wg_ind <dbl>, ws_ind <dbl>, 조절변인 <lgl>,
## # N3_3_15_b_분류 <dbl>, 가구소득1M <dbl>, 가구소득2M <dbl>,
## # 가구소득3M <dbl>, 가구소득4M <dbl>, 종속변인 <lgl>, SGDS_K총점 <dbl>,
## # 독립변인 <lgl>, 영양총점 <dbl>
sdata <- rdata[,c(29:34,42:45,47,49)]
head(sdata)
## # A tibble: 6 x 12
## ID_A3 ID_A4_2 ID_A5 ID_A6_1 ID_A7 B1 가구소득1M 가구소득2M 가구소득3M
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 70 4 5 1 3 0 1 0
## 2 1 71 2 7 2 2 -1 -1 -1
## 3 2 71 2 5 2 2 -1 -1 -1
## 4 1 82 2 7 2 4 0 0 0
## 5 1 68 2 3 1 2 -1 -1 -1
## 6 2 66 2 3 1 2 -1 -1 -1
## # ... with 3 more variables: 가구소득4M <dbl>, SGDS_K총점 <dbl>,
## # 영양총점 <dbl>
colnames(sdata) <- c("c1","c2","c3","c4","c5","c6","m1","m2","m3","m4","y","x")
head(sdata)
## # A tibble: 6 x 12
## c1 c2 c3 c4 c5 c6 m1 m2 m3 m4 y x
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 70 4 5 1 3 0 1 0 0 10 14
## 2 1 71 2 7 2 2 -1 -1 -1 -1 2 5
## 3 2 71 2 5 2 2 -1 -1 -1 -1 2 3
## 4 1 82 2 7 2 4 0 0 0 1 7 5
## 5 1 68 2 3 1 2 -1 -1 -1 -1 3 2
## 6 2 66 2 3 1 2 -1 -1 -1 -1 2 2
#covariate
model1 <- lm(y ~ c1+c2+c3+c4+c5+c6, data = sdata)
summary(model1)
##
## Call:
## lm(formula = y ~ c1 + c2 + c3 + c4 + c5 + c6, data = sdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8838 -2.4293 -0.6742 2.0948 13.6242
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.728269 0.496525 -11.537 < 2e-16 ***
## c1 -0.326929 0.080046 -4.084 4.46e-05 ***
## c2 0.048603 0.006011 8.086 6.88e-16 ***
## c3 0.458176 0.061685 7.428 1.19e-13 ***
## c4 -0.359390 0.027279 -13.175 < 2e-16 ***
## c5 0.909174 0.079011 11.507 < 2e-16 ***
## c6 1.738008 0.037602 46.221 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.497 on 10076 degrees of freedom
## (216 observations deleted due to missingness)
## Multiple R-squared: 0.2729, Adjusted R-squared: 0.2725
## F-statistic: 630.3 on 6 and 10076 DF, p-value: < 2.2e-16
#x + covariate
model2 <- lm(y ~ c1+c2+c3+c4+c5+c6+x, data = sdata)
summary(model2)
##
## Call:
## lm(formula = y ~ c1 + c2 + c3 + c4 + c5 + c6 + x, data = sdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8238 -2.3127 -0.5646 1.9843 13.3239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.940139 0.477835 -8.246 < 2e-16 ***
## c1 -0.176135 0.076625 -2.299 0.0215 *
## c2 0.029035 0.005777 5.026 5.09e-07 ***
## c3 0.048182 0.060390 0.798 0.4250
## c4 -0.253026 0.026285 -9.626 < 2e-16 ***
## c5 0.813202 0.075547 10.764 < 2e-16 ***
## c6 1.319390 0.038368 34.388 < 2e-16 ***
## x 0.418797 0.013482 31.062 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.341 on 10075 degrees of freedom
## (216 observations deleted due to missingness)
## Multiple R-squared: 0.3364, Adjusted R-squared: 0.336
## F-statistic: 729.8 on 7 and 10075 DF, p-value: < 2.2e-16
#x + covariate + moderation
model3 <- lm(y ~ c1+c2+c3+c4+c5+c6+x+m1+m2+m3+m4, data = sdata)
summary(model3)
##
## Call:
## lm(formula = y ~ c1 + c2 + c3 + c4 + c5 + c6 + x + m1 + m2 +
## m3 + m4, data = sdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.7030 -2.3067 -0.5607 1.9589 13.3762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.5670385 0.4848724 -7.357 2.03e-13 ***
## c1 -0.1801004 0.0766066 -2.351 0.018743 *
## c2 0.0256534 0.0058182 4.409 1.05e-05 ***
## c3 0.0003381 0.0621251 0.005 0.995658
## c4 -0.2341856 0.0267264 -8.762 < 2e-16 ***
## c5 0.7928578 0.0756893 10.475 < 2e-16 ***
## c6 1.3183782 0.0383431 34.384 < 2e-16 ***
## x 0.4080190 0.0136840 29.817 < 2e-16 ***
## m1 0.2600182 0.0708848 3.668 0.000246 ***
## m2 0.1340920 0.0657536 2.039 0.041445 *
## m3 -0.0500651 0.0663185 -0.755 0.450315
## m4 -0.1839538 0.0692013 -2.658 0.007867 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.338 on 10071 degrees of freedom
## (216 observations deleted due to missingness)
## Multiple R-squared: 0.3379, Adjusted R-squared: 0.3372
## F-statistic: 467.2 on 11 and 10071 DF, p-value: < 2.2e-16
#x + covariate + moderation + interaction
model4 <- lm(y ~ c1+c2+c3+c4+c5+c6+x+m1+m2+m3+m4+x*m1+x*m2+x*m3+x*m4, data = sdata)
summary(model4)
##
## Call:
## lm(formula = y ~ c1 + c2 + c3 + c4 + c5 + c6 + x + m1 + m2 +
## m3 + m4 + x * m1 + x * m2 + x * m3 + x * m4, data = sdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.7912 -2.2920 -0.5648 1.9529 13.3842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.615575 0.485039 -7.454 9.79e-14 ***
## c1 -0.182620 0.076608 -2.384 0.0172 *
## c2 0.026307 0.005822 4.518 6.31e-06 ***
## c3 0.001468 0.062123 0.024 0.9811
## c4 -0.237275 0.026826 -8.845 < 2e-16 ***
## c5 0.798335 0.075703 10.546 < 2e-16 ***
## c6 1.323863 0.038406 34.470 < 2e-16 ***
## x 0.396771 0.014568 27.236 < 2e-16 ***
## m1 0.024641 0.118960 0.207 0.8359
## m2 0.091791 0.110651 0.830 0.4068
## m3 0.036732 0.107417 0.342 0.7324
## m4 -0.016443 0.108764 -0.151 0.8798
## x:m1 0.053534 0.021350 2.507 0.0122 *
## x:m2 0.017983 0.023912 0.752 0.4520
## x:m3 -0.019000 0.025481 -0.746 0.4559
## x:m4 -0.050568 0.027912 -1.812 0.0701 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.337 on 10067 degrees of freedom
## (216 observations deleted due to missingness)
## Multiple R-squared: 0.3385, Adjusted R-squared: 0.3375
## F-statistic: 343.4 on 15 and 10067 DF, p-value: < 2.2e-16
6. Interaction
library(effects)
## Warning: package 'effects' was built under R version 3.5.3
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
interm1 <- effect("x*m1", model4, xlevels=list(x=c(0,5,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
inter1 <- as.data.frame(interm1)
head(inter1)
## x m1 fit se lower upper
## 1 0 -1 2.654757 0.13558114 2.388991 2.920523
## 2 5 -1 4.367811 0.09168271 4.188095 4.547527
## 3 10 -1 6.080865 0.20267201 5.683587 6.478142
## 4 0 0 2.679398 0.05955740 2.562654 2.796142
## 5 5 0 4.660121 0.04235924 4.577088 4.743153
## 6 10 0 6.640843 0.10272796 6.439476 6.842210
interm2 <- effect("x*m2", model4, xlevels=list(x=c(0,5,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
inter2 <- as.data.frame(interm2)
head(inter2)
## x m2 fit se lower upper
## 1 0 -1.0 2.585243 0.12882983 2.332711 2.837776
## 2 5 -1.0 4.485014 0.09147725 4.305700 4.664327
## 3 10 -1.0 6.384784 0.21002662 5.973090 6.796478
## 4 0 -0.5 2.631139 0.08376633 2.466940 2.795338
## 5 5 -0.5 4.575867 0.06098163 4.456331 4.695403
## 6 10 -0.5 6.520595 0.14282817 6.240624 6.800567
interm3 <- effect("x*m3", model4, xlevels=list(x=c(0,5,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
inter3 <- as.data.frame(interm3)
head(inter3)
## x m3 fit se lower upper
## 1 0 -1.0 2.642587 0.12643346 2.394752 2.890422
## 2 5 -1.0 4.733860 0.09431136 4.548991 4.918729
## 3 10 -1.0 6.825134 0.21662270 6.400510 7.249758
## 4 0 -0.5 2.660953 0.08300204 2.498252 2.823653
## 5 5 -0.5 4.704728 0.06105492 4.585048 4.824407
## 6 10 -0.5 6.748502 0.14176723 6.470610 7.026394
interm4 <- effect("x*m4", model4, xlevels=list(x=c(0,5,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
inter4 <- as.data.frame(interm4)
head(inter4)
## x m4 fit se lower upper
## 1 0 -1.0 2.697204 0.12641844 2.449399 2.945010
## 2 5 -1.0 4.946575 0.09749707 4.755461 5.137688
## 3 10 -1.0 7.195945 0.22344314 6.757952 7.633938
## 4 0 -0.5 2.688983 0.08251835 2.527230 2.850735
## 5 5 -0.5 4.811933 0.06016340 4.694001 4.929865
## 6 10 -0.5 6.934883 0.13851257 6.663371 7.206395
inter1$m1 <- factor(inter1$m1, levels = c(-1, 0, 1), labels = c("Grand mean", "Othes quintile", "Quintile 1"))
inter1$x <- factor(inter1$x, levels = c(0, 5, 10), labels = c("Good", "warning", "Bad"))
inter2$m2 <- factor(inter2$m2, levels = c(-1, 0, 1), labels = c("Grand mean", "Othes quintile", "Quintile 2"))
inter2$x <- factor(inter2$x, levels = c(0, 5, 10), labels = c("Good", "warning", "Bad"))
inter3$m3 <- factor(inter3$m3, levels = c(-1, 0, 1), labels = c("Grand mean", "Othes quintile", "Quintile 3"))
inter3$x <- factor(inter3$x, levels = c(0, 5, 10), labels = c("Good", "warning", "Bad"))
inter4$m4 <- factor(inter4$m4, levels = c(-1, 0, 1), labels = c("Grand mean", "Othes quintile", "Quintile 4"))
inter4$x <- factor(inter4$x, levels = c(0, 5, 10), labels = c("Good", "warning", "Bad"))
7. Interaction plot
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
plot.int1 <- ggplot(data = inter1, aes(x = x, y = fit, group = m1)) +
geom_line(size = 2, aes(color = m1)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 1")
rinter2 <- inter2[-c(4:6,10:12),]
plot.int2 <- ggplot(data = rinter2, aes(x = x, y = fit, group = m2)) +
geom_line(size = 2, aes(color = m2)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 2")
rinter3 <- inter3[-c(4:6,10:12),]
plot.int3 <- ggplot(data = rinter3, aes(x = x, y = fit, group = m3)) +
geom_line(size = 2, aes(color = m3)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 3")
rinter4 <- inter4[-c(4:6,10:12),]
plot.int4 <- ggplot(data = rinter4, aes(x = x, y = fit, group = m4)) +
geom_line(size = 2, aes(color = m4)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 4")
library(ggpubr)
ggarrange(plot.int1,plot.int2,plot.int3,plot.int4)
plot.int1
plot.int2
plot.int3
plot.int4
8. New Interaction plot 1
#Interaction
library(effects)
ninterm1 <- effect("x*m1", model4, xlevels=list(x=c(0,10), m1=c(-1,1)), se = T, confidence.level = .95, typical = mean)
ninter1 <- as.data.frame(ninterm1)
head(ninter1)
## x m1 fit se lower upper
## 1 0 -1 2.654757 0.1355811 2.388991 2.920523
## 2 10 -1 6.080865 0.2026720 5.683587 6.478142
## 3 0 1 2.704039 0.1304414 2.448348 2.959730
## 4 10 1 7.200821 0.1349983 6.936198 7.465445
ninterm2 <- effect("x*m2", model4, xlevels=list(x=c(0,10), m1=c(-1,1)), se = T, confidence.level = .95, typical = mean)
ninter2 <- as.data.frame(ninterm2)
head(ninter2)
## x m2 fit se lower upper
## 1 0 -1.0 2.585243 0.12882983 2.332711 2.837776
## 2 10 -1.0 6.384784 0.21002662 5.973090 6.796478
## 3 0 -0.5 2.631139 0.08376633 2.466940 2.795338
## 4 10 -0.5 6.520595 0.14282817 6.240624 6.800567
## 5 0 0.0 2.677034 0.05965141 2.560106 2.793963
## 6 10 0.0 6.656407 0.10138608 6.457670 6.855144
ninterm3 <- effect("x*m3", model4, xlevels=list(x=c(0,10), m1=c(-1,1)), se = T, confidence.level = .95, typical = mean)
ninter3 <- as.data.frame(ninterm3)
head(ninter3)
## x m3 fit se lower upper
## 1 0 -1.0 2.642587 0.12643346 2.394752 2.890422
## 2 10 -1.0 6.825134 0.21662270 6.400510 7.249758
## 3 0 -0.5 2.660953 0.08300204 2.498252 2.823653
## 4 10 -0.5 6.748502 0.14176723 6.470610 7.026394
## 5 0 0.0 2.679319 0.05968612 2.562322 2.796316
## 6 10 0.0 6.671871 0.10064701 6.474582 6.869159
ninterm4 <- effect("x*m4", model4, xlevels=list(x=c(0,10), m1=c(-1,1)), se = T, confidence.level = .95, typical = mean)
ninter4 <- as.data.frame(ninterm4)
head(ninter4)
## x m4 fit se lower upper
## 1 0 -1.0 2.697204 0.12641844 2.449399 2.945010
## 2 10 -1.0 7.195945 0.22344314 6.757952 7.633938
## 3 0 -0.5 2.688983 0.08251835 2.527230 2.850735
## 4 10 -0.5 6.934883 0.13851257 6.663371 7.206395
## 5 0 0.0 2.680761 0.05959633 2.563940 2.797582
## 6 10 0.0 6.673821 0.10006735 6.477669 6.869973
ninter1$m1 <- factor(ninter1$m1, levels = c(-1, 1), labels = c("Grand mean", "Quintile 1"))
ninter1$x <- factor(ninter1$x, levels = c(0, 10), labels = c("Good", "Bad"))
ninter2$m2 <- factor(ninter2$m2, levels = c(-1, 1), labels = c("Grand mean", "Quintile 2"))
ninter2$x <- factor(ninter2$x, levels = c(0, 10), labels = c("Good", "Bad"))
ninter3$m3 <- factor(ninter3$m3, levels = c(-1, 1), labels = c("Grand mean", "Quintile 3"))
ninter3$x <- factor(ninter3$x, levels = c(0, 10), labels = c("Good", "Bad"))
ninter4$m4 <- factor(ninter4$m4, levels = c(-1, 1), labels = c("Grand mean", "Quintile 4"))
ninter4$x <- factor(ninter4$x, levels = c(0, 10), labels = c("Good", "Bad"))
#Interaction plot
library(ggplot2)
nplot.int1 <- ggplot(data = ninter1, aes(x = x, y = fit, group = m1)) +
geom_line(size = 2, aes(color = m1)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 1")
nrinter2 <- ninter2[-c(3:8),]
nplot.int2 <- ggplot(data = nrinter2, aes(x = x, y = fit, group = m2)) +
geom_line(size = 2, aes(color = m2)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 2")
nrinter3 <- ninter3[-c(3:8),]
nplot.int3 <- ggplot(data = nrinter3, aes(x = x, y = fit, group = m3)) +
geom_line(size = 2, aes(color = m3)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 3")
nrinter4 <- ninter4[-c(3:8),]
nplot.int4 <- ggplot(data = nrinter4, aes(x = x, y = fit, group = m4)) +
geom_line(size = 2, aes(color = m4)) +
ylim(2, 8) +
ylab("Depression Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 4")
library(ggpubr)
ggarrange(nplot.int1,nplot.int2,nplot.int3,nplot.int4)
nplot.int1
nplot.int2
nplot.int3
nplot.int4
9. New Interaction plot 2
library(effects)
nninterm1 <- effect("x*m1", model4, xlevels=list(x=c(0,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
nninter1 <- as.data.frame(nninterm1)
head(nninter1)
## x m1 fit se lower upper
## 1 0 -1 2.654757 0.1355811 2.388991 2.920523
## 2 10 -1 6.080865 0.2026720 5.683587 6.478142
## 3 0 0 2.679398 0.0595574 2.562654 2.796142
## 4 10 0 6.640843 0.1027280 6.439476 6.842210
## 5 0 1 2.704039 0.1304414 2.448348 2.959730
## 6 10 1 7.200821 0.1349983 6.936198 7.465445
nninterm2 <- effect("x*m2", model4, xlevels=list(x=c(0,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
nninter2 <- as.data.frame(nninterm2)
head(nninter2)
## x m2 fit se lower upper
## 1 0 -1.0 2.585243 0.12882983 2.332711 2.837776
## 2 10 -1.0 6.384784 0.21002662 5.973090 6.796478
## 3 0 -0.5 2.631139 0.08376633 2.466940 2.795338
## 4 10 -0.5 6.520595 0.14282817 6.240624 6.800567
## 5 0 0.0 2.677034 0.05965141 2.560106 2.793963
## 6 10 0.0 6.656407 0.10138608 6.457670 6.855144
nninterm3 <- effect("x*m3", model4, xlevels=list(x=c(0,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
nninter3 <- as.data.frame(nninterm3)
head(nninter3)
## x m3 fit se lower upper
## 1 0 -1.0 2.642587 0.12643346 2.394752 2.890422
## 2 10 -1.0 6.825134 0.21662270 6.400510 7.249758
## 3 0 -0.5 2.660953 0.08300204 2.498252 2.823653
## 4 10 -0.5 6.748502 0.14176723 6.470610 7.026394
## 5 0 0.0 2.679319 0.05968612 2.562322 2.796316
## 6 10 0.0 6.671871 0.10064701 6.474582 6.869159
nninterm4 <- effect("x*m4", model4, xlevels=list(x=c(0,10), m1=c(-1,0,1)), se = T, confidence.level = .95, typical = mean)
nninter4 <- as.data.frame(nninterm4)
head(nninter4)
## x m4 fit se lower upper
## 1 0 -1.0 2.697204 0.12641844 2.449399 2.945010
## 2 10 -1.0 7.195945 0.22344314 6.757952 7.633938
## 3 0 -0.5 2.688983 0.08251835 2.527230 2.850735
## 4 10 -0.5 6.934883 0.13851257 6.663371 7.206395
## 5 0 0.0 2.680761 0.05959633 2.563940 2.797582
## 6 10 0.0 6.673821 0.10006735 6.477669 6.869973
nninter1$m1 <- factor(nninter1$m1, levels = c(-1, 0, 1), labels = c("Grand mean", "Others quintile", "Quintile 1"))
nninter1$x <- factor(nninter1$x, levels = c(0, 10), labels = c("Good", "Bad"))
nninter2$m2 <- factor(nninter2$m2, levels = c(-1, 0, 1), labels = c("Grand mean", "Others quintile", "Quintile 2"))
nninter2$x <- factor(nninter2$x, levels = c(0, 10), labels = c("Good", "Bad"))
nninter3$m3 <- factor(nninter3$m3, levels = c(-1, 0, 1), labels = c("Grand mean", "Others quintile", "Quintile 3"))
nninter3$x <- factor(nninter3$x, levels = c(0, 10), labels = c("Good", "Bad"))
nninter4$m4 <- factor(nninter4$m4, levels = c(-1, 0, 1), labels = c("Grand mean", "Others quintile", "Quintile 4"))
nninter4$x <- factor(nninter4$x, levels = c(0, 10), labels = c("Good", "Bad"))
#Interaction plot
library(ggplot2)
nnplot.int1 <- ggplot(data = nninter1, aes(x = x, y = fit, group = m1)) +
geom_line(size = 1, aes(color = m1)) +
ylim(2, 8) +
ylab("Depressive Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 1")
nnrinter2 <- nninter2[-c(3:4,7:8),]
nnplot.int2 <- ggplot(data = nnrinter2, aes(x = x, y = fit, group = m2)) +
geom_line(size = 1, aes(color = m2)) +
ylim(2, 8) +
ylab("Depressive Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 2")
nnrinter3 <- nninter3[-c(3:4,7:8),]
nnplot.int3 <- ggplot(data = nnrinter3, aes(x = x, y = fit, group = m3)) +
geom_line(size = 1, aes(color = m3)) +
ylim(2, 8) +
ylab("Depressive Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 3")
nnrinter4 <- nninter4[-c(3:4,7:8),]
nnplot.int4 <- ggplot(data = nnrinter4, aes(x = x, y = fit, group = m4)) +
geom_line(size = 1, aes(color = m4)) +
ylim(2, 8) +
ylab("Depressive Symptoms") +
xlab("Nutrition management status") +
ggtitle("Quintile 4")
library(ggpubr)
ggarrange(nnplot.int1,nnplot.int2,nnplot.int3,nnplot.int4)
nnplot.int1
nnplot.int2
nnplot.int3
nnplot.int4
10. Dummy variable
newdata <- rdata[,c(29:34,41,47,49)]
head(newdata)
## # A tibble: 6 x 9
## ID_A3 ID_A4_2 ID_A5 ID_A6_1 ID_A7 B1 N3_3_15_b_분류 SGDS_K총점
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 70 4 5 1 3 2 10
## 2 1 71 2 7 2 2 5 2
## 3 2 71 2 5 2 2 5 2
## 4 1 82 2 7 2 4 4 7
## 5 1 68 2 3 1 2 5 3
## 6 2 66 2 3 1 2 5 2
## # ... with 1 more variable: 영양총점 <dbl>
colnames(newdata) <- c("c1","c2","c3","c4","c5","c6","income","y","x")
head(newdata)
## # A tibble: 6 x 9
## c1 c2 c3 c4 c5 c6 income y x
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 70 4 5 1 3 2 10 14
## 2 1 71 2 7 2 2 5 2 5
## 3 2 71 2 5 2 2 5 2 3
## 4 1 82 2 7 2 4 4 7 5
## 5 1 68 2 3 1 2 5 3 2
## 6 2 66 2 3 1 2 5 2 2
ddata <- newdata %>% mutate(d1 = ifelse(income == 5, 0,
ifelse(income == 4, 1, 0)),
d2 = ifelse(income == 5, 0,
ifelse(income == 3, 1, 0)),
d3 = ifelse(income == 5, 0,
ifelse(income == 2, 1, 0)),
d4 = ifelse(income == 5, 0,
ifelse(income == 1, 1, 0)))
dummyfit <- lm(y ~ c1+c2+c3+c4+c5+c6+x+x*d1+x*d2+x*d3+x*d4, data = ddata)
summary(dummyfit)
##
## Call:
## lm(formula = y ~ c1 + c2 + c3 + c4 + c5 + c6 + x + x * d1 + x *
## d2 + x * d3 + x * d4, data = ddata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.7912 -2.2920 -0.5648 1.9529 13.3842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.752296 0.494244 -7.592 3.43e-14 ***
## c1 -0.182620 0.076608 -2.384 0.0172 *
## c2 0.026307 0.005822 4.518 6.31e-06 ***
## c3 0.001468 0.062123 0.024 0.9811
## c4 -0.237275 0.026826 -8.845 < 2e-16 ***
## c5 0.798335 0.075703 10.546 < 2e-16 ***
## c6 1.323863 0.038406 34.470 < 2e-16 ***
## x 0.394822 0.036896 10.701 < 2e-16 ***
## d1 0.120277 0.170546 0.705 0.4807
## d2 0.173453 0.170676 1.016 0.3095
## d3 0.228512 0.175200 1.304 0.1922
## d4 0.161362 0.184348 0.875 0.3814
## x:d1 -0.048619 0.047915 -1.015 0.3103
## x:d2 -0.017050 0.045739 -0.373 0.7093
## x:d3 0.019933 0.044327 0.450 0.6530
## x:d4 0.055483 0.042104 1.318 0.1876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.337 on 10067 degrees of freedom
## (216 observations deleted due to missingness)
## Multiple R-squared: 0.3385, Adjusted R-squared: 0.3375
## F-statistic: 343.4 on 15 and 10067 DF, p-value: < 2.2e-16