library(readxl)
library(dplyr)
library(magrittr)
library(knitr)
library(writexl)
library(agricolae)
library(DescTools)
library(dunn.test)read_excel("~/Desktop/raw_winter.xlsx", sheet = "raw") -> raw_winter
read_excel("~/Desktop/raw_winter.xlsx", sheet = "KRDI") -> ppKRDI
read_excel("~/Desktop/raw_winter.xlsx", sheet = "FRAS4") -> ppFRAS4
raw_winter %>%
dplyr::mutate(KRDA = kcal / RDA) -> raw_winter
raw_winter$agesex <- factor(raw_winter$agesex, levels=c("M_Young", "F_Young", "M_Old", "F_Old"))
raw_winter$prepost <- factor(raw_winter$prepost, levels=c("pre", "post"))
raw_winter## # A tibble: 222 × 19
## ...1 ...2 agesex offin prepost sex age carrer height weight carbohydr…¹
## <dbl> <chr> <fct> <chr> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2014 F_Old off <NA> F 21 12 158 50.6 96.4
## 2 2 2014 M_Old off <NA> M 21 11 168. 62.4 186.
## 3 3 2014 M_Old off <NA> M 20 14 180. 78.3 62.7
## 4 4 2014 F_Old off <NA> F 19 10 163. 60.0 310.
## 5 5 2014 M_Old off <NA> M 23 14 178. 74.8 177.
## 6 6 2014 M_Old off <NA> M 20 11 184. 75.4 252.
## 7 7 2014 M_Old off <NA> M 19 8 180. 76.5 209.
## 8 9 2014 M_Old off <NA> M 20 12 184. 77 181.
## 9 10 2014 M_Old off <NA> M 19 11 168 63.9 171.
## 10 11 2014 M_Old off <NA> M 21 15 172. 63.8 192.
## # … with 212 more rows, 8 more variables: `fat(g)` <dbl>, `protein(g)` <dbl>,
## # `dietary fibre(g)` <dbl>, kcal <dbl>, RDA <dbl>, droms <dbl>, BAP <dbl>,
## # KRDA <dbl>, and abbreviated variable name ¹`carbohydrate(g)`raw_winter %>%
group_by(agesex) %>%
summarise(n = n(),
Age = mean(age),
Age_sd = sd(age),
Weight= mean(weight),
Weight_sd = sd(weight),
Krda = mean(KRDA)*100,
Krda_sd = sd(KRDA)*100) -> a
raw_winter %>%
tidyr::drop_na(carrer) -> carrerfull
carrerfull %>%
group_by(agesex) %>%
summarise(n = n(),
Carrer = mean(carrer),
Carrer_sd = sd(carrer)) -> b
raw_winter %>%
tidyr::drop_na(height) -> heightfull
heightfull %>%
group_by(agesex) %>%
summarise(n = n(),
Height = mean(height),
Height_sd = sd(height)) -> c
df1 <- data.frame(c(a,b,c))
df1## agesex n Age Age_sd Weight Weight_sd Krda Krda_sd agesex.1
## 1 M_Young 99 15.19192 2.023720 59.05909 13.768602 88.20138 40.68680 M_Young
## 2 F_Young 60 14.96667 1.775222 53.56333 7.262732 92.26304 34.34027 F_Young
## 3 M_Old 45 21.60000 2.434599 68.75889 5.422829 76.10686 30.73211 M_Old
## 4 F_Old 18 21.22222 2.734433 59.36111 4.285425 68.89143 50.68055 F_Old
## n.1 Carrer Carrer_sd agesex.2 n.2 Height Height_sd
## 1 67 6.447761 2.996380 M_Young 68 167.4294 8.018766
## 2 42 6.515952 2.909439 F_Young 42 159.8548 6.240005
## 3 38 13.000000 3.653950 M_Old 38 174.3000 5.674552
## 4 14 12.285714 3.383866 F_Old 14 164.0000 3.268968raw_winter %>%
tidyr::drop_na(offin) -> offinfull
offinfull$offin <- factor(offinfull$offin, levels=c("off", "in"))
offinfull %>%
group_by(offin) %>%
summarise(n = n(),
Krda = mean(KRDA)*100,
Krda_sd = sd(KRDA)*100) -> d
df2 <- data.frame(d)
df2## offin n Krda Krda_sd
## 1 off 139 89.80567 42.02902
## 2 in 58 75.93017 29.28151raw_winter %>%
group_by(prepost) %>%
summarise(n = n(),
Age = mean(age),
Age_sd = sd(age),
Carrer = mean(carrer),
Carrer_sd = sd(carrer),
Height = mean(height),
Height_sd = sd(height),
Weight= mean(weight),
Weight_sd = sd(weight),
Krda = mean(KRDA)*100,
Krda_sd = sd(KRDA)*100)-> e
ppFRAS4 %>%
tidyr::drop_na(c(pre.droms, pre.BAP, post.droms, post.BAP)) -> ppFRAS4
ppFRAS4 %>%
summarise(n = n(),
pre.d.roms = mean(pre.droms),
pre.d.roms_sd = sd(pre.droms),
post.d.roms = mean(post.droms),
post.d.roms_sd = sd(post.droms),
pre.Bap = mean(pre.BAP),
pre.Bap_sd = sd(pre.BAP),
post.Bap = mean(post.BAP),
post.Bap_sd = sd(post.BAP)) -> f
df3 <- as.data.frame(e)
df3## prepost n Age Age_sd Carrer Carrer_sd Height Height_sd Weight
## 1 pre 11 17.36364 4.924890 9.090909 4.504543 166.5727 7.967947 58.05455
## 2 post 11 18.36364 4.924890 10.090909 4.504543 167.7273 8.580104 58.15455
## 3 <NA> 200 16.81500 3.404221 NA NA NA NA 59.72500
## Weight_sd Krda Krda_sd
## 1 9.495406 87.73899 28.62289
## 2 9.189272 76.77473 23.71104
## 3 11.753193 85.61461 39.76818df4 <- as.data.frame(f)
df4## n pre.d.roms pre.d.roms_sd post.d.roms post.d.roms_sd pre.Bap pre.Bap_sd
## 1 9 236.3333 37.71604 198.3333 41.88675 2623.778 497.9171
## post.Bap post.Bap_sd
## 1 2456.778 728.7981# writexl::write_xlsx(df,"~/Desktop/result_winter.xlsx")
# writexl::write_xlsx(setNames(df1,"setting"), paste0(outpath, "test.xlsx"))
# writexl::write_xlsx(list("temperature" = df1, "humidity"=df2), paste0(outpath, "test.xlsx"))
write_xlsx(list("DS_nutrition" = df1, "DS_offin"=df2, "DS_prepost"=df3, "FRAS4" =df4), paste0("~/Desktop/result_winter.xlsx"))by(raw_winter$KRDA, raw_winter$agesex, shapiro.test)## raw_winter$agesex: M_Young
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.65674, p-value = 7.005e-14
##
## ------------------------------------------------------------
## raw_winter$agesex: F_Young
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96607, p-value = 0.09346
##
## ------------------------------------------------------------
## raw_winter$agesex: M_Old
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.93161, p-value = 0.01069
##
## ------------------------------------------------------------
## raw_winter$agesex: F_Old
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.82722, p-value = 0.00378M_Young : p-value = 7.005e-14 < 유의수준 = 0.05, 정규성 가정을
만족하지 않음
F_Young : p-value = 0.09346 > 유의수준 = 0.05, 정규성 가정을
만족함
M_Old : p-value = 0.01069 < 유의수준 = 0.05, 정규성 가정을
만족하지 않음
F_Old : p-value = 0.00378 < 유의수준 = 0.05, 정규성 가정을
만족하지 않음
정규성 가정을 전부 다 만족하지 않으므로 비모수검정
kruskal.test(formula = KRDA ~ agesex,
data = raw_winter)##
## Kruskal-Wallis rank sum test
##
## data: KRDA by agesex
## Kruskal-Wallis chi-squared = 13.654, df = 3, p-value = 0.003417# summary(aov(KRDA ~ agesex, raw_winter))
# anova(lm(KRDA ~ agesex, raw_winter))p-value = 0.003417 < 유의수준 = 0.05, 통계적으로 유의한 차이 증명.
## bonferroni -------------------------------------------------------
dunn.test(raw_winter$KRDA, raw_winter$agesex, method = 'bonferroni')## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 13.6537, df = 3, p-value = 0
##
##
## Comparison of x by group
## (Bonferroni)
## Col Mean-|
## Row Mean | F_Old F_Young M_Old
## ---------+---------------------------------
## F_Young | -3.098451
## | 0.0058*
## |
## M_Old | -1.177301 2.557522
## | 0.7172 0.0316
## |
## M_Young | -2.648306 0.941847 -1.948176
## | 0.0243* 1.0000 0.1542
##
## alpha = 0.05
## Reject Ho if p <= alpha/2## Tukey ------------------------------------------------------------
nparcomp::nparcomp(formula = KRDA ~ agesex,
data = raw_winter,
type = "Tukey")##
## #------Nonparametric Multiple Comparisons for relative contrast effects-----#
##
## - Alternative Hypothesis: True relative contrast effect p is not equal to 1/2
## - Type of Contrast : Tukey
## - Confidence level: 95 %
## - Method = Logit - Transformation
## - Estimation Method: Pairwise rankings
##
## #---------------------------Interpretation----------------------------------#
## p(a,b) > 1/2 : b tends to be larger than a
## #---------------------------------------------------------------------------#
## ## $Data.Info
## Sample Size
## 1 M_Young 99
## 2 F_Young 60
## 3 M_Old 45
## 4 F_Old 18
##
## $Contrast
## M_Young F_Young M_Old F_Old
## F_Young - M_Young -1 1 0 0
## M_Old - M_Young -1 0 1 0
## F_Old - M_Young -1 0 0 1
## M_Old - F_Young 0 -1 1 0
## F_Old - F_Young 0 -1 0 1
## F_Old - M_Old 0 0 -1 1
##
## $Analysis
## Comparison Estimator Lower Upper Statistic p.Value
## 1 p( M_Young , F_Young ) 0.553 0.428 0.671 1.070137 0.69697249
## 2 p( M_Young , M_Old ) 0.389 0.262 0.533 -1.957642 0.18579773
## 3 p( M_Young , F_Old ) 0.300 0.128 0.556 -1.992360 0.17335663
## 4 p( F_Young , M_Old ) 0.360 0.236 0.507 -2.409155 0.06644819
## 5 p( F_Young , F_Old ) 0.289 0.126 0.534 -2.197727 0.11120924
## 6 p( M_Old , F_Old ) 0.373 0.179 0.619 -1.308648 0.53923776
##
## $Overall
## Quantile p.Value
## 1 2.525806 0.06644819
##
## $input
## $input$formula
## KRDA ~ agesex
##
## $input$data
## # A tibble: 222 × 19
## ...1 ...2 agesex offin prepost sex age carrer height weight carbohydr…¹
## <dbl> <chr> <fct> <chr> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2014 F_Old off <NA> F 21 12 158 50.6 96.4
## 2 2 2014 M_Old off <NA> M 21 11 168. 62.4 186.
## 3 3 2014 M_Old off <NA> M 20 14 180. 78.3 62.7
## 4 4 2014 F_Old off <NA> F 19 10 163. 60.0 310.
## 5 5 2014 M_Old off <NA> M 23 14 178. 74.8 177.
## 6 6 2014 M_Old off <NA> M 20 11 184. 75.4 252.
## 7 7 2014 M_Old off <NA> M 19 8 180. 76.5 209.
## 8 9 2014 M_Old off <NA> M 20 12 184. 77 181.
## 9 10 2014 M_Old off <NA> M 19 11 168 63.9 171.
## 10 11 2014 M_Old off <NA> M 21 15 172. 63.8 192.
## # … with 212 more rows, 8 more variables: `fat(g)` <dbl>, `protein(g)` <dbl>,
## # `dietary fibre(g)` <dbl>, kcal <dbl>, RDA <dbl>, droms <dbl>, BAP <dbl>,
## # KRDA <dbl>, and abbreviated variable name ¹`carbohydrate(g)`
##
## $input$type
## [1] "Tukey"
##
## $input$conf.level
## [1] 0.95
##
## $input$alternative
## [1] "two.sided" "less" "greater"
##
## $input$asy.method
## [1] "logit" "probit" "normal" "mult.t"
##
## $input$plot.simci
## [1] FALSE
##
## $input$control
## NULL
##
## $input$info
## [1] TRUE
##
## $input$rounds
## [1] 3
##
## $input$contrast.matrix
## NULL
##
## $input$correlation
## [1] FALSE
##
## $input$weight.matrix
## [1] FALSE
##
##
## $text.Output
## [1] "True relative contrast effect p is not equal to 1/2"
##
## $connames
## [1] "p( M_Young , F_Young )" "p( M_Young , M_Old )" "p( M_Young , F_Old )"
## [4] "p( F_Young , M_Old )" "p( F_Young , F_Old )" "p( M_Old , F_Old )"
##
## $AsyMethod
## [1] "Logit - Transformation"
##
## attr(,"class")
## [1] "nparcomp"## LSD --------------------------------------------------------------
model<-aov(KRDA ~ agesex, raw_winter)
comparison<-LSD.test(model,"agesex",group=F)
comparison## $statistics
## MSerror Df Mean CV
## 0.1454256 218 0.8528186 44.71609
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none agesex 4 0.05
##
## $means
## KRDA std r LCL UCL Min Max Q25
## F_Old 0.6889143 0.5068055 18 0.5117606 0.8660679 0.1455912 2.225234 0.3444273
## F_Young 0.9226304 0.3434027 60 0.8255994 1.0196614 0.3626002 1.987853 0.6537550
## M_Old 0.7610686 0.3073211 45 0.6490268 0.8731104 0.3088165 1.753853 0.5327090
## M_Young 0.8820138 0.4068680 99 0.8064753 0.9575524 0.3268118 3.993003 0.6812488
## Q50 Q75
## F_Old 0.4958059 0.946188
## F_Young 0.8960822 1.140565
## M_Old 0.6594340 1.002902
## M_Young 0.8108828 1.021395
##
## $comparison
## difference pvalue signif. LCL UCL
## F_Old - F_Young -0.23371614 0.0235 * -0.43570235 -0.031729939
## F_Old - M_Old -0.07215432 0.4982 -0.28176531 0.137456669
## F_Old - M_Young -0.19309957 0.0494 * -0.38568586 -0.000513286
## F_Young - M_Old 0.16156182 0.0328 * 0.01334447 0.309779174
## F_Young - M_Young 0.04061657 0.5157 -0.08235129 0.163584430
## M_Old - M_Young -0.12094525 0.0791 . -0.25607273 0.014182225
##
## $groups
## NULL
##
## attr(,"class")
## [1] "group"comparison<-LSD.test(model,"agesex",p.adj="bonferroni",group=F)
comparison## $statistics
## MSerror Df Mean CV
## 0.1454256 218 0.8528186 44.71609
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD bonferroni agesex 4 0.05
##
## $means
## KRDA std r LCL UCL Min Max Q25
## F_Old 0.6889143 0.5068055 18 0.5117606 0.8660679 0.1455912 2.225234 0.3444273
## F_Young 0.9226304 0.3434027 60 0.8255994 1.0196614 0.3626002 1.987853 0.6537550
## M_Old 0.7610686 0.3073211 45 0.6490268 0.8731104 0.3088165 1.753853 0.5327090
## M_Young 0.8820138 0.4068680 99 0.8064753 0.9575524 0.3268118 3.993003 0.6812488
## Q50 Q75
## F_Old 0.4958059 0.946188
## F_Young 0.8960822 1.140565
## M_Old 0.6594340 1.002902
## M_Young 0.8108828 1.021395
##
## $comparison
## difference pvalue signif. LCL UCL
## F_Old - F_Young -0.23371614 0.1413 -0.50658472 0.03915244
## F_Old - M_Old -0.07215432 1.0000 -0.35532343 0.21101479
## F_Old - M_Young -0.19309957 0.2964 -0.45326956 0.06707041
## F_Young - M_Old 0.16156182 0.1967 -0.03866897 0.36179262
## F_Young - M_Young 0.04061657 1.0000 -0.12550401 0.20673715
## M_Old - M_Young -0.12094525 0.4747 -0.30349259 0.06160208
##
## $groups
## NULL
##
## attr(,"class")
## [1] "group"## scheffe ; group 별 n수 다를 때 -----------------------------------
scheffe.test(model, "agesex", alpha = 0.05, console = T)##
## Study: model ~ "agesex"
##
## Scheffe Test for KRDA
##
## Mean Square Error : 0.1454256
##
## agesex, means
##
## KRDA std r Min Max
## F_Old 0.6889143 0.5068055 18 0.1455912 2.225234
## F_Young 0.9226304 0.3434027 60 0.3626002 1.987853
## M_Old 0.7610686 0.3073211 45 0.3088165 1.753853
## M_Young 0.8820138 0.4068680 99 0.3268118 3.993003
##
## Alpha: 0.05 ; DF Error: 218
## Critical Value of F: 2.646014
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Means with the same letter are not significantly different.
##
## KRDA groups
## F_Young 0.9226304 a
## M_Young 0.8820138 a
## M_Old 0.7610686 a
## F_Old 0.6889143 a## bonferroni -------------------------------------------------------
PostHocTest(model, method='bonferroni')##
## Posthoc multiple comparisons of means : Bonferroni
## 95% family-wise confidence level
##
## $agesex
## diff lwr.ci upr.ci pval
## F_Young-M_Young 0.04061657 -0.1255040 0.20673715 1.0000
## M_Old-M_Young -0.12094525 -0.3034926 0.06160208 0.4747
## F_Old-M_Young -0.19309957 -0.4532696 0.06707041 0.2964
## M_Old-F_Young -0.16156182 -0.3617926 0.03866897 0.1967
## F_Old-F_Young -0.23371614 -0.5065847 0.03915244 0.1413
## F_Old-M_Old -0.07215432 -0.3553234 0.21101479 1.0000
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1by(raw_winter$KRDA, raw_winter$offin, shapiro.test)## raw_winter$offin: in
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96033, p-value = 0.05542
##
## ------------------------------------------------------------
## raw_winter$offin: off
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.78062, p-value = 3.799e-13in : p-value = 0.05542 > 유의수준 = 0.05, 정규성 가정을
만족함
off : p-value = 3.799e-13 < 유의수준 = 0.05, 정규성 가정을
만족하지 않음
정규성 가정을 전부 다 만족하지 않으므로 비모수검정
wilcox.test(formula = KRDA ~ offin,
data = raw_winter,
alternative = "two.sided")##
## Wilcoxon rank sum test with continuity correction
##
## data: KRDA by offin
## W = 3125, p-value = 0.01304
## alternative hypothesis: true location shift is not equal to 0p-value = 0.01304 < 유의수준 = 0.05, 통계적으로 유의한 차이 증명.
by(raw_winter$KRDA, raw_winter$prepost, shapiro.test)## raw_winter$prepost: pre
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.98767, p-value = 0.9941
##
## ------------------------------------------------------------
## raw_winter$prepost: post
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.97595, p-value = 0.9395pre : p-value = 0.9941 > 유의수준 = 0.05, 정규성 가정을
만족함
post : p-value = 0.9395 > 유의수준 = 0.05, 정규성 가정을
만족함
정규성 가정을 전부 다 만족하므로 모수검정
by(raw_winter$droms, raw_winter$prepost, shapiro.test)## raw_winter$prepost: pre
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.93018, p-value = 0.4497
##
## ------------------------------------------------------------
## raw_winter$prepost: post
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.93455, p-value = 0.4941pre : p-value = 0.4497 > 유의수준 = 0.05, 정규성 가정을
만족함
post : p-value = 0.4941 > 유의수준 = 0.05, 정규성 가정을
만족함
정규성 가정을 전부 다 만족하므로 모수검정
by(raw_winter$BAP, raw_winter$prepost, shapiro.test)## raw_winter$prepost: pre
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.75673, p-value = 0.00431
##
## ------------------------------------------------------------
## raw_winter$prepost: post
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.80918, p-value = 0.01243pre : p-value = 0.00431 < 유의수준 = 0.05, 정규성 가정을
만족하지 않음
post : p-value = 0.01243 < 유의수준 = 0.05, 정규성 가정을
만족하지 않음
정규성 가정을 전부 다 만족하지 않으므로 비모수검정
#KRDA
t.test(ppKRDI$pre,
ppKRDI$post,
alternative = "two.sided",
paired = TRUE)##
## Paired t-test
##
## data: ppKRDI$pre and ppKRDI$post
## t = 1.9068, df = 10, p-value = 0.08565
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.01847447 0.23775952
## sample estimates:
## mean difference
## 0.1096425#d-roms
t.test(ppFRAS4$pre.droms,
ppFRAS4$post.droms,
alternative = "two.sided",
paired = TRUE)##
## Paired t-test
##
## data: ppFRAS4$pre.droms and ppFRAS4$post.droms
## t = 4.2007, df = 8, p-value = 0.002994
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 17.13948 58.86052
## sample estimates:
## mean difference
## 38#BAP
wilcox.test(ppFRAS4$pre.BAP,
ppFRAS4$post.BAP,
alternative = "two.sided",
paired = FALSE)## Warning in wilcox.test.default(ppFRAS4$pre.BAP, ppFRAS4$post.BAP, alternative =
## "two.sided", : cannot compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: ppFRAS4$pre.BAP and ppFRAS4$post.BAP
## W = 44.5, p-value = 0.7572
## alternative hypothesis: true location shift is not equal to 0KRDA : p-value = 0.08565 > 유의수준 = 0.05, 통계적으로 유의하지 않음.
d-roms : p-value = 0.002994 < 유의수준 = 0.05, 통계적으로 유의한 차이 증명.
BAP : p-value = 0.7572 > 유의수준 = 0.05, 통계적으로 유의하지 않음.