## >> CARREGAR PACOTES
library(tidyverse)
## -- Attaching packages --------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0 v purrr 0.2.5
## v tibble 1.4.2 v dplyr 0.7.8
## v tidyr 0.8.1 v stringr 1.3.1
## v readr 1.1.1 v forcats 0.3.0
## -- Conflicts ------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:purrr':
##
## transpose
library(corrplot) # https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
## corrplot 0.84 loaded
library(agricolae)
library(NCStats) # para instalar, rode o comando: source("http://www.rforge.net/NCStats/InstallNCStats.R")
## Loading required package: FSA
## ## FSA v0.8.22. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.
## ## NCStats v0.4.7 by Derek H. Ogle, Northland College.
##
## Attaching package: 'NCStats'
## The following object is masked from 'package:FSA':
##
## rSquared
library(qwraps2)
library(summarytools)
##
## Attaching package: 'summarytools'
## The following object is masked from 'package:NCStats':
##
## view
## >> CARREGAR BANCO DE DADOS
banco = readstata13::read.dta13(file.choose()) # bancofinalwide281118.dta
## >> RENOMEAR VARIÁVEIS
banco=rename(banco,
parity=nativivoss1,
schooling=anos_de_estudo1,
mother_age=idade_mae_calc1,
total_energy=new_kcal_total2,
total_protein=prottotal2,
total_fat=liptotal2,
total_carbohydrates=carbtotal2,
total_fiber=fibrtotal2,
saturated_fat=sattotal2,
unsaturated_fat=instotal2,
calcium=calciototal2,
iron=ferrototal2,
tryptophan=triptofanototal2,
HMOSecretor=HMOSecretor3,
Diversity_7d=Diversity3,
Evenness_7d=Evenness3,
"2'FL_7d"=FL2nmolmL3,
"3FL_7d"=FL3nmolmL3,
LNnT_7d=LNnTnmolmL3,
"3'SL_7d"=SL3nmolmL3,
DFLac_7d=DFLacnmolmL3,
"6'SL_7d"=SL6nmolmL3,
LNT_7d=LNTnmolmL3,
LNFPI_7d=LNFPInmolmL3,
LNFPII_7d=LNFPIInmolmL3,
LNFPIII_7d=LNFPIIInmolmL3,
LSTb_7d=LSTbnmolmL3,
LSTc_7d=LSTcnmolmL3,
DFLNT_7d=DFLNTnmolmL3,
LNH_7d=LNHnmolmL3,
DSLNT_7d=DSLNTnmolmL3,
FLNH_7d=FLNHnmolmL3,
DFLNH_7d=DFLNHnmolmL3,
FDSLNH_7d=FDSLNHnmolmL3,
DSLNH_7d=DSLNHnmolmL3,
SUM_7d=SUMnmolmL3,
Sia_7d=SianmolmL3,
Fuc_7d=FucnmolmL3,
HMOSecretor_1m=HMOSecretor4,
Diversity_1m=Diversity4,
Evenness_1m=Evenness4,
"2'FL_1m"=FL2nmolmL4,
"3FL_1m"=FL3nmolmL4,
LNnT_1m=LNnTnmolmL4,
"3'SL_1m"=SL3nmolmL4,
DFLac_1m=DFLacnmolmL4,
"6'SL_1m"=SL6nmolmL4,
LNT_1m=LNTnmolmL4,
LNFPI_1m=LNFPInmolmL4,
LNFPII_1m=LNFPIInmolmL4,
LNFPIII_1m=LNFPIIInmolmL4,
LSTb_1m=LSTbnmolmL4,
LSTc_1m=LSTcnmolmL4,
DFLNT_1m=DFLNTnmolmL4,
LNH_1m=LNHnmolmL4,
DSLNT_1m=DSLNTnmolmL4,
FLNH_1m=FLNHnmolmL4,
DFLNH_1m=DFLNHnmolmL4,
FDSLNH_1m=FDSLNHnmolmL4,
DSLNH_1m=DSLNHnmolmL4,
SUM_1m=SUMnmolmL4,
Sia_1m=SianmolmL4,
Fuc_1m=FucnmolmL4,
HMOSecretor_3m=HMOSecretor5,
Diversity_3m=Diversity5,
Evenness_3m=Evenness5,
"2'FL_3m"=FL2nmolmL5,
"3FL_3m"=FL3nmolmL5,
LNnT_3m=LNnTnmolmL5,
"3'SL_3m"=SL3nmolmL5,
DFLac_3m=DFLacnmolmL5,
"6'SL_3m"=SL6nmolmL5,
LNT_3m=LNTnmolmL5,
LNFPI_3m=LNFPInmolmL5,
LNFPII_3m=LNFPIInmolmL5,
LNFPIII_3m=LNFPIIInmolmL5,
LSTb_3m=LSTbnmolmL5,
LSTc_3m=LSTcnmolmL5,
DFLNT_3m=DFLNTnmolmL5,
LNH_3m=LNHnmolmL5,
DSLNT_3m=DSLNTnmolmL5,
FLNH_3m=FLNHnmolmL5,
DFLNH_3m=DFLNHnmolmL5,
FDSLNH_3m=FDSLNHnmolmL5,
DSLNH_3m=DSLNHnmolmL5,
SUM_3m=SUMnmolmL5,
Sia_3m=SianmolmL5,
Fuc_3m=FucnmolmL5)
##############################################
#### 1-FACTOR ANOVA + Bonferroni #############
##############################################
## PREPARAR...
tb2 = banco %>% select("record_id", 238:257, 583:602, 928:947) #dim = 148 x 60 (=20 x 3)
tb2 = melt(tb2, id="record_id")
headtail(tb2)
## record_id variable value
## 1 1 2'FL_7d NA
## 2 2 2'FL_7d NA
## 3 3 2'FL_7d NA
## 8878 146 SUM_3m NA
## 8879 147 SUM_3m NA
## 8880 148 SUM_3m NA
tb2$onda=rep(c("w1","w2","w3"), each=2960)
tb2$onda=as.factor(tb2$onda)
headtail(tb2)
## record_id variable value onda
## 1 1 2'FL_7d NA w1
## 2 2 2'FL_7d NA w1
## 3 3 2'FL_7d NA w1
## 8878 146 SUM_3m NA w3
## 8879 147 SUM_3m NA w3
## 8880 148 SUM_3m NA w3
tb2$variable = as.factor(sub("[_].*", "", tb2$variable)) # levels(tb2$variable)
headtail(tb2)
## record_id variable value onda
## 1 1 2'FL NA w1
## 2 2 2'FL NA w1
## 3 3 2'FL NA w1
## 8878 146 SUM NA w3
## 8879 147 SUM NA w3
## 8880 148 SUM NA w3
## APONTAR...
mean_sd_w1 = tapply(tb2$value[tb2$onda=="w1"], tb2$variable[tb2$onda=="w1"], mean_sd, na_rm=T, show_n="never")
mean_sd_w2 = tapply(tb2$value[tb2$onda=="w2"], tb2$variable[tb2$onda=="w2"], mean_sd, na_rm=T, show_n="never")
mean_sd_w3 = tapply(tb2$value[tb2$onda=="w3"], tb2$variable[tb2$onda=="w3"], mean_sd, na_rm=T, show_n="never")
knitr::kable(mean_sd_w1, format="pandoc", col.names=c("Mean (SD)"), caption = "Title of the table")
Title of the table
| 2’FL |
4,575.30 \(\pm\) 2,394.27 |
| 3’SL |
342.46 \(\pm\) 187.23 |
| 3FL |
235.14 \(\pm\) 116.77 |
| 6’SL |
350.21 \(\pm\) 224.24 |
| DFLac |
253.30 \(\pm\) 134.79 |
| DFLNH |
160.79 \(\pm\) 147.04 |
| DFLNT |
1,451.19 \(\pm\) 790.97 |
| DSLNH |
237.30 \(\pm\) 104.67 |
| DSLNT |
635.73 \(\pm\) 218.85 |
| FDSLNH |
110.29 \(\pm\) 102.04 |
| FLNH |
119.35 \(\pm\) 82.04 |
| LNFPI |
2,638.96 \(\pm\) 1,377.29 |
| LNFPII |
1,047.32 \(\pm\) 653.09 |
| LNFPIII |
99.20 \(\pm\) 50.33 |
| LNH |
88.89 \(\pm\) 93.16 |
| LNnT |
814.53 \(\pm\) 486.69 |
| LNT |
1,617.76 \(\pm\) 892.83 |
| LSTb |
120.78 \(\pm\) 108.50 |
| LSTc |
825.28 \(\pm\) 321.24 |
| SUM |
15,723.77 \(\pm\) 2,636.81 |
knitr::kable(mean_sd_w2, format="markdown", col.names=c("Mean (SD)"))
| 2’FL |
5,001.34 \(\pm\) 2,487.45 |
| 3’SL |
440.48 \(\pm\) 219.58 |
| 3FL |
335.93 \(\pm\) 174.35 |
| 6’SL |
578.32 \(\pm\) 227.09 |
| DFLac |
385.91 \(\pm\) 186.37 |
| DFLNH |
172.78 \(\pm\) 124.98 |
| DFLNT |
1,434.74 \(\pm\) 836.71 |
| DSLNH |
287.09 \(\pm\) 123.18 |
| DSLNT |
264.65 \(\pm\) 221.52 |
| FDSLNH |
283.94 \(\pm\) 220.80 |
| FLNH |
135.71 \(\pm\) 86.44 |
| LNFPI |
1,721.06 \(\pm\) 1,233.38 |
| LNFPII |
1,308.19 \(\pm\) 721.88 |
| LNFPIII |
77.10 \(\pm\) 39.09 |
| LNH |
108.42 \(\pm\) 78.93 |
| LNnT |
352.32 \(\pm\) 224.29 |
| LNT |
1,496.85 \(\pm\) 650.19 |
| LSTb |
123.09 \(\pm\) 49.69 |
| LSTc |
295.55 \(\pm\) 128.66 |
| SUM |
14,803.45 \(\pm\) 2,402.59 |
knitr::kable(mean_sd_w3, format="markdown", col.names=c("Mean (SD)"))
| 2’FL |
5,937.33 \(\pm\) 3,257.14 |
| 3’SL |
672.92 \(\pm\) 434.00 |
| 3FL |
373.56 \(\pm\) 155.32 |
| 6’SL |
344.50 \(\pm\) 131.90 |
| DFLac |
425.27 \(\pm\) 177.84 |
| DFLNH |
65.18 \(\pm\) 55.60 |
| DFLNT |
1,327.04 \(\pm\) 718.03 |
| DSLNH |
128.27 \(\pm\) 54.49 |
| DSLNT |
57.79 \(\pm\) 74.31 |
| FDSLNH |
355.06 \(\pm\) 259.80 |
| FLNH |
88.59 \(\pm\) 52.32 |
| LNFPI |
1,164.30 \(\pm\) 838.89 |
| LNFPII |
1,745.22 \(\pm\) 808.22 |
| LNFPIII |
52.00 \(\pm\) 24.67 |
| LNH |
70.60 \(\pm\) 40.95 |
| LNnT |
672.56 \(\pm\) 287.32 |
| LNT |
1,412.98 \(\pm\) 674.99 |
| LSTb |
95.00 \(\pm\) 38.70 |
| LSTc |
131.58 \(\pm\) 78.46 |
| SUM |
15,119.75 \(\pm\) 2,554.96 |
#by(tb2$value[tb2$onda=="w1"], tb2$variable[tb2$onda=="w1"], mean_sd, na_rm=T, show_n="never")
## FOGOOO !!!
levels(tb2$variable)
## [1] "2'FL" "3'SL" "3FL" "6'SL" "DFLac" "DFLNH" "DFLNT"
## [8] "DSLNH" "DSLNT" "FDSLNH" "FLNH" "LNFPI" "LNFPII" "LNFPIII"
## [15] "LNH" "LNnT" "LNT" "LSTb" "LSTc" "SUM"
model_1=aov(value[variable=="2'FL"] ~ onda[variable=="2'FL"], data=tb2)
out_1=LSD.test(model_1, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_1 ~ "onda"
##
## LSD t Test for value[variable == "2'FL"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 6785050
##
## onda[variable == "2'FL"], means and individual ( 95 %) CI
##
## value.variable.....2.FL.. std r LCL UCL Min
## w1 4575.304 2394.269 43 3788.539 5362.068 8.873210
## w2 5001.336 2487.450 55 4305.675 5696.997 1.473807
## w3 5937.334 3257.144 21 4811.513 7063.155 47.255995
## Max
## w1 9436.263
## w2 10619.831
## w3 12839.160
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -426.0320 1.0000 -1714.092 862.0281
## w1 - w3 -1362.0298 0.1557 -3046.583 322.5237
## w2 - w3 -935.9978 0.4918 -2559.132 687.1367
#out_1
model_2=aov(value[variable=="3'SL"] ~ onda[variable=="3'SL"], data=tb2)
out_2=LSD.test(model_2, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_2 ~ "onda"
##
## LSD t Test for value[variable == "3'SL"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 67611.51
##
## onda[variable == "3'SL"], means and individual ( 95 %) CI
##
## value.variable.....3.SL.. std r LCL UCL Min
## w1 342.4587 187.2320 43 263.9210 420.9965 49.48814
## w2 440.4754 219.5762 55 371.0319 509.9189 141.66936
## w3 672.9190 433.9960 21 560.5353 785.3026 192.10786
## Max
## w1 965.1158
## w2 1014.8000
## w3 1771.0216
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -98.01669 0.1998 -226.5957 30.56227
## w1 - w3 -330.46024 0.0000 *** -498.6186 -162.30185
## w2 - w3 -232.44355 0.0021 ** -394.4709 -70.41623
#out_2
model_3=aov(value[variable=="3FL"] ~ onda[variable=="3FL"], data=tb2)
out_3=LSD.test(model_3, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_3 ~ "onda"
##
## LSD t Test for value[variable == "3FL"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 23246.63
##
## onda[variable == "3FL"], means and individual ( 95 %) CI
##
## value.variable.....3FL.. std r LCL UCL Min
## w1 235.1399 116.7747 43 189.0879 281.1918 15.47513
## w2 335.9340 174.3458 55 295.2147 376.6534 102.52644
## w3 373.5552 155.3181 21 307.6572 439.4533 145.13431
## Max
## w1 531.6593
## w2 1090.3054
## w3 708.7105
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -100.7942 0.0046 ** -176.1887 -25.39970
## w1 - w3 -138.4154 0.0027 ** -237.0179 -39.81282
## w2 - w3 -37.6212 1.0000 -132.6287 57.38630
#out_3
model_4=aov(value[variable=="6'SL"] ~ onda[variable=="6'SL"], data=tb2)
out_4=LSD.test(model_4, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_4 ~ "onda"
##
## LSD t Test for value[variable == "6'SL"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 45213.68
##
## onda[variable == "6'SL"], means and individual ( 95 %) CI
##
## value.variable.....6.SL.. std r LCL UCL Min
## w1 350.2109 224.2430 43 285.9861 414.4358 66.30064
## w2 578.3153 227.0937 55 521.5274 635.1033 193.59157
## w3 344.5038 131.9007 21 252.6012 436.4064 100.23206
## Max
## w1 1253.9314
## w2 1168.9649
## w3 532.4109
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -228.104400 0e+00 *** -333.2508 -122.9580
## w1 - w3 5.707099 1e+00 -131.8057 143.2199
## w2 - w3 233.811499 1e-04 *** 101.3124 366.3106
#out_4
model_5=aov(value[variable=="DFLac"] ~ onda[variable=="DFLac"], data=tb2)
out_5=LSD.test(model_5, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_5 ~ "onda"
##
## LSD t Test for value[variable == "DFLac"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 28200.08
##
## onda[variable == "DFLac"], means and individual ( 95 %) CI
##
## value.variable.....DFLac.. std r LCL UCL Min
## w1 253.2968 134.7861 43 202.5751 304.0184 8.383564
## w2 385.9064 186.3712 55 341.0580 430.7547 24.976855
## w3 425.2714 177.8392 21 352.6912 497.8515 10.772385
## Max
## w1 601.4064
## w2 872.8305
## w3 759.9184
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -132.6096 5e-04 *** -215.6491 -49.57005
## w1 - w3 -171.9746 6e-04 *** -280.5755 -63.37366
## w2 - w3 -39.3650 1e+00 -144.0063 65.27630
#out_5
model_6=aov(value[variable=="DFLNH"] ~ onda[variable=="DFLNH"], data=tb2)
out_6=LSD.test(model_6, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_6 ~ "onda"
##
## LSD t Test for value[variable == "DFLNH"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 15632.47
##
## onda[variable == "DFLNH"], means and individual ( 95 %) CI
##
## value.variable.....DFLNH.. std r LCL UCL Min
## w1 160.78542 147.03719 43 123.02107 198.5498 7.736311
## w2 172.77919 124.98161 55 139.38776 206.1706 7.324023
## w3 65.17786 55.60009 21 11.13895 119.2168 17.595097
## Max
## w1 793.9381
## w2 606.3908
## w3 206.9968
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -11.99377 1.0000 -73.82008 49.83255
## w1 - w3 95.60756 0.0145 * 14.74974 176.46538
## w2 - w3 107.60132 0.0032 ** 29.69159 185.51106
#out_6
model_7=aov(value[variable=="DFLNT"] ~ onda[variable=="DFLNT"], data=tb2)
out_7=LSD.test(model_7, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_7 ~ "onda"
##
## LSD t Test for value[variable == "DFLNT"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 641314.1
##
## onda[variable == "DFLNT"], means and individual ( 95 %) CI
##
## value.variable.....DFLNT.. std r LCL UCL Min
## w1 1451.191 790.9748 43 1209.3087 1693.073 148.66005
## w2 1434.740 836.7051 55 1220.8670 1648.614 13.26986
## w3 1327.039 718.0333 21 980.9174 1673.160 64.93601
## Max
## w1 3532.190
## w2 3379.887
## w3 2270.841
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 16.45048 1 -379.5494 412.4503
## w1 - w3 124.15215 1 -393.7452 642.0495
## w2 - w3 107.70168 1 -391.3130 606.7164
#out_7
model_8=aov(value[variable=="DSLNH"] ~ onda[variable=="DSLNH"], data=tb2)
out_8=LSD.test(model_8, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_8 ~ "onda"
##
## LSD t Test for value[variable == "DSLNH"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 11541.81
##
## onda[variable == "DSLNH"], means and individual ( 95 %) CI
##
## value.variable.....DSLNH.. std r LCL UCL Min
## w1 237.2962 104.67281 43 204.84685 269.7454 57.39176
## w2 287.0922 123.17568 55 258.40041 315.7841 89.36123
## w3 128.2675 54.48875 21 81.83412 174.7008 14.77203
## Max
## w1 479.5164
## w2 556.1867
## w3 247.4448
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -49.79608 0.0739 . -102.92080 3.328639
## w1 - w3 109.02869 0.0007 *** 39.55101 178.506373
## w2 - w3 158.82477 0.0000 *** 91.88026 225.769290
#out_8
model_9=aov(value[variable=="DSLNT"] ~ onda[variable=="DSLNT"], data=tb2)
out_9=LSD.test(model_9, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_9 ~ "onda"
##
## LSD t Test for value[variable == "DSLNT"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 41138.35
##
## onda[variable == "DSLNT"], means and individual ( 95 %) CI
##
## value.variable.....DSLNT.. std r LCL UCL Min
## w1 635.72806 218.85386 43 574.46599 696.9901 25.979887
## w2 264.64741 221.52380 55 210.47919 318.8156 3.194394
## w3 57.79312 74.31118 21 -29.86987 145.4561 7.785786
## Max
## w1 1130.3700
## w2 736.7102
## w3 336.1953
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 371.0806 0e+00 *** 270.78478 471.3765
## w1 - w3 577.9349 0e+00 *** 446.76579 709.1041
## w2 - w3 206.8543 4e-04 *** 80.46759 333.2410
#out_9
model_10=aov(value[variable=="FDSLNH"] ~ onda[variable=="FDSLNH"], data=tb2)
out_10=LSD.test(model_10, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_10 ~ "onda"
##
## LSD t Test for value[variable == "FDSLNH"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 38102.92
##
## onda[variable == "FDSLNH"], means and individual ( 95 %) CI
##
## value.variable.....FDSLNH.. std r LCL UCL Min
## w1 110.2931 102.0421 43 51.33442 169.2517 20.18971
## w2 283.9391 220.8017 55 231.80755 336.0706 37.33255
## w3 355.0608 259.8006 21 270.69389 439.4276 87.96471
## Max
## w1 478.1588
## w2 992.3354
## w3 1092.6768
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -173.6460 0.0001 *** -270.1708 -77.12124
## w1 - w3 -244.7677 0.0000 *** -371.0049 -118.53050
## w2 - w3 -71.1217 0.4745 -192.7563 50.51288
#out_10
model_11=aov(value[variable=="FLNH"] ~ onda[variable=="FLNH"], data=tb2)
out_11=LSD.test(model_11, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_11 ~ "onda"
##
## LSD t Test for value[variable == "FLNH"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 6387.157
##
## onda[variable == "FLNH"], means and individual ( 95 %) CI
##
## value.variable.....FLNH.. std r LCL UCL Min
## w1 119.3451 82.03619 43 95.20598 143.4843 10.23279
## w2 135.7129 86.44340 55 114.36893 157.0568 15.16136
## w3 88.5876 52.31639 21 54.04566 123.1295 17.50512
## Max
## w1 363.1713
## w2 408.1923
## w3 237.2466
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -16.36775 0.9494 -55.887434 23.15193
## w1 - w3 30.75754 0.4529 -20.927170 82.44225
## w2 - w3 47.12529 0.0699 . -2.674985 96.92557
#out_11
model_12=aov(value[variable=="LNFPI"] ~ onda[variable=="LNFPI"], data=tb2)
out_12=LSD.test(model_12, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_12 ~ "onda"
##
## LSD t Test for value[variable == "LNFPI"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 1516313
##
## onda[variable == "LNFPI"], means and individual ( 95 %) CI
##
## value.variable.....LNFPI.. std r LCL UCL Min
## w1 2638.964 1377.2946 43 2267.0327 3010.895 57.73936
## w2 1721.061 1233.3797 55 1392.1979 2049.925 92.17351
## w3 1164.302 838.8868 21 632.0865 1696.517 125.01014
## Max
## w1 5000.461
## w2 5171.549
## w3 2965.465
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 917.9027 0.0011 ** 308.9914 1526.814
## w1 - w3 1474.6623 0.0000 *** 678.3147 2271.010
## w2 - w3 556.7596 0.2418 -210.5530 1324.072
#out_12
model_13=aov(value[variable=="LNFPII"] ~ onda[variable=="LNFPII"], data=tb2)
out_13=LSD.test(model_13, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_13 ~ "onda"
##
## LSD t Test for value[variable == "LNFPII"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 509643.2
##
## onda[variable == "LNFPII"], means and individual ( 95 %) CI
##
## value.variable.....LNFPII.. std r LCL UCL Min
## w1 1047.319 653.0938 43 831.6926 1262.945 267.5823
## w2 1308.185 721.8783 55 1117.5278 1498.843 328.8927
## w3 1745.219 808.2214 21 1436.6687 2053.769 735.7669
## Max
## w1 2921.199
## w2 3433.966
## w3 3627.576
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -260.8667 0.2257 -613.8813 92.147940
## w1 - w3 -697.9002 0.0011 ** -1159.5805 -236.219908
## w2 - w3 -437.0335 0.0559 . -881.8809 7.813869
#out_13
model_14=aov(value[variable=="LNFPIII"] ~ onda[variable=="LNFPIII"], data=tb2)
out_14=LSD.test(model_14, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_14 ~ "onda"
##
## LSD t Test for value[variable == "LNFPIII"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 1733.484
##
## onda[variable == "LNFPIII"], means and individual ( 95 %) CI
##
## value.variable.....LNFPIII.. std r LCL UCL Min
## w1 99.20448 50.32939 43 86.62889 111.78006 45.24402
## w2 77.09704 39.09277 55 65.97765 88.21643 23.76174
## w3 51.99749 24.66874 21 34.00245 69.99253 17.37602
## Max
## w1 300.21614
## w2 248.89678
## w3 89.47754
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 22.10744 0.0309 * 1.5191818 42.69569
## w1 - w3 47.20699 0.0001 *** 20.2812139 74.13276
## w2 - w3 25.09955 0.0614 . -0.8445058 51.04361
#out_14
model_15=aov(value[variable=="LNH"] ~ onda[variable=="LNH"], data=tb2)
out_15=LSD.test(model_15, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_15 ~ "onda"
##
## LSD t Test for value[variable == "LNH"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 6331.485
##
## onda[variable == "LNH"], means and individual ( 95 %) CI
##
## value.variable.....LNH.. std r LCL UCL Min
## w1 88.89172 93.16362 43 64.85800 112.9254 9.598484
## w2 108.42429 78.92617 55 87.17356 129.6750 21.820551
## w3 70.59994 40.94547 21 36.20887 104.9910 16.490096
## Max
## w1 525.6633
## w2 428.2573
## w3 175.0604
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -19.53257 0.6909 -58.87964 19.81450
## w1 - w3 18.29178 1.0000 -33.16718 69.75075
## w2 - w3 37.82435 0.1992 -11.75841 87.40712
#out_15
model_16=aov(value[variable=="LNnT"] ~ onda[variable=="LNnT"], data=tb2)
out_16=LSD.test(model_16, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_16 ~ "onda"
##
## LSD t Test for value[variable == "LNnT"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 123413.3
##
## onda[variable == "LNnT"], means and individual ( 95 %) CI
##
## value.variable.....LNnT.. std r LCL UCL Min
## w1 814.5270 486.6904 43 708.4187 920.6353 154.00172
## w2 352.3151 224.2884 55 258.4937 446.1366 67.55773
## w3 672.5642 287.3165 21 520.7286 824.3999 324.57999
## Max
## w1 2048.859
## w2 1118.779
## w3 1498.809
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 462.2119 0.0000 *** 288.49560 635.9282
## w1 - w3 141.9628 0.3953 -85.22723 369.1528
## w2 - w3 -320.2491 0.0017 ** -539.15574 -101.3425
#out_16
model_17=aov(value[variable=="LNT"] ~ onda[variable=="LNT"], data=tb2)
out_17=LSD.test(model_17, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_17 ~ "onda"
##
## LSD t Test for value[variable == "LNT"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 563970.9
##
## onda[variable == "LNT"], means and individual ( 95 %) CI
##
## value.variable.....LNT.. std r LCL UCL Min
## w1 1617.760 892.8264 43 1390.932 1844.588 286.4739
## w2 1496.848 650.1947 55 1296.286 1697.410 542.2827
## w3 1412.980 674.9859 21 1088.401 1737.560 519.0140
## Max
## w1 5211.305
## w2 3721.284
## w3 3167.202
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 120.91225 1.0000 -250.4416 492.2661
## w1 - w3 204.77988 0.9235 -280.8849 690.4447
## w2 - w3 83.86763 1.0000 -384.0898 551.8250
#out_17
model_18=aov(value[variable=="LSTb"] ~ onda[variable=="LSTb"], data=tb2)
out_18=LSD.test(model_18, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_18 ~ "onda"
##
## LSD t Test for value[variable == "LSTb"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 5669.823
##
## onda[variable == "LSTb"], means and individual ( 95 %) CI
##
## value.variable.....LSTb.. std r LCL UCL Min
## w1 120.77566 108.49924 43 98.03239 143.5189 37.49274
## w2 123.09114 49.68634 55 102.98143 143.2009 38.22470
## w3 94.99957 38.70427 21 62.45507 127.5441 31.99769
## Max
## w1 760.9990
## w2 308.9165
## w3 167.3872
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 -2.315477 1.0000 -39.54988 34.91892
## w1 - w3 25.776090 0.6032 -22.91988 74.47206
## w2 - w3 28.091567 0.4456 -18.82894 75.01208
#out_18
model_19=aov(value[variable=="LSTc"] ~ onda[variable=="LSTc"], data=tb2)
out_19=LSD.test(model_19, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_19 ~ "onda"
##
## LSD t Test for value[variable == "LSTc"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 46129.97
##
## onda[variable == "LSTc"], means and individual ( 95 %) CI
##
## value.variable.....LSTc.. std r LCL UCL Min
## w1 825.2800 321.23687 43 760.40764 890.1524 201.83792
## w2 295.5469 128.65637 55 238.18646 352.9074 79.63516
## w3 131.5765 78.46418 21 38.74736 224.4057 23.57499
## Max
## w1 1692.6517
## w2 668.2915
## w3 384.6444
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 529.7331 0.0000 *** 423.52655 635.9396
## w1 - w3 693.7035 0.0000 *** 554.80426 832.6027
## w2 - w3 163.9704 0.0107 * 30.13546 297.8054
#out_19
model_20=aov(value[variable=="SUM"] ~ onda[variable=="SUM"], data=tb2)
out_20=LSD.test(model_20, "onda", p.adj="bonferroni", console=T, group=F)
##
## Study: model_20 ~ "onda"
##
## LSD t Test for value[variable == "SUM"]
## P value adjustment method: bonferroni
##
## Mean Square Error: 6330049
##
## onda[variable == "SUM"], means and individual ( 95 %) CI
##
## value.variable.....SUM.. std r LCL UCL Min
## w1 15723.77 2636.812 43 14963.84 16483.70 9099.194
## w2 14803.45 2402.593 55 14131.52 15475.38 8241.586
## w3 15119.75 2554.962 21 14032.33 16207.16 8583.907
## Max
## w1 18503.51
## w2 18243.79
## w3 19008.19
##
## Alpha: 0.05 ; DF Error: 116
## Critical Value of t: 2.429194
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## w1 - w2 920.3232 0.2248 -323.7992 2164.446
## w1 - w3 604.0234 1.0000 -1023.0674 2231.114
## w2 - w3 -316.2999 1.0000 -1884.0667 1251.467
#out_20
print(out_1$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -426.0320 1.0000 -1714.092 862.0281
## w1 - w3 -1362.0298 0.1557 -3046.583 322.5237
## w2 - w3 -935.9978 0.4918 -2559.132 687.1367
print(out_2$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -98.01669 0.1998 -226.5957 30.56227
## w1 - w3 -330.46024 0.0000 *** -498.6186 -162.30185
## w2 - w3 -232.44355 0.0021 ** -394.4709 -70.41623
print(out_3$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -100.7942 0.0046 ** -176.1887 -25.39970
## w1 - w3 -138.4154 0.0027 ** -237.0179 -39.81282
## w2 - w3 -37.6212 1.0000 -132.6287 57.38630
print(out_4$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -228.104400 0e+00 *** -333.2508 -122.9580
## w1 - w3 5.707099 1e+00 -131.8057 143.2199
## w2 - w3 233.811499 1e-04 *** 101.3124 366.3106
print(out_5$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -132.6096 5e-04 *** -215.6491 -49.57005
## w1 - w3 -171.9746 6e-04 *** -280.5755 -63.37366
## w2 - w3 -39.3650 1e+00 -144.0063 65.27630
print(out_6$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -11.99377 1.0000 -73.82008 49.83255
## w1 - w3 95.60756 0.0145 * 14.74974 176.46538
## w2 - w3 107.60132 0.0032 ** 29.69159 185.51106
print(out_7$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 16.45048 1 -379.5494 412.4503
## w1 - w3 124.15215 1 -393.7452 642.0495
## w2 - w3 107.70168 1 -391.3130 606.7164
print(out_8$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -49.79608 0.0739 . -102.92080 3.328639
## w1 - w3 109.02869 0.0007 *** 39.55101 178.506373
## w2 - w3 158.82477 0.0000 *** 91.88026 225.769290
print(out_9$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 371.0806 0e+00 *** 270.78478 471.3765
## w1 - w3 577.9349 0e+00 *** 446.76579 709.1041
## w2 - w3 206.8543 4e-04 *** 80.46759 333.2410
print(out_10$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -173.6460 0.0001 *** -270.1708 -77.12124
## w1 - w3 -244.7677 0.0000 *** -371.0049 -118.53050
## w2 - w3 -71.1217 0.4745 -192.7563 50.51288
print(out_11$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -16.36775 0.9494 -55.887434 23.15193
## w1 - w3 30.75754 0.4529 -20.927170 82.44225
## w2 - w3 47.12529 0.0699 . -2.674985 96.92557
print(out_12$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 917.9027 0.0011 ** 308.9914 1526.814
## w1 - w3 1474.6623 0.0000 *** 678.3147 2271.010
## w2 - w3 556.7596 0.2418 -210.5530 1324.072
print(out_13$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -260.8667 0.2257 -613.8813 92.147940
## w1 - w3 -697.9002 0.0011 ** -1159.5805 -236.219908
## w2 - w3 -437.0335 0.0559 . -881.8809 7.813869
print(out_14$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 22.10744 0.0309 * 1.5191818 42.69569
## w1 - w3 47.20699 0.0001 *** 20.2812139 74.13276
## w2 - w3 25.09955 0.0614 . -0.8445058 51.04361
print(out_15$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -19.53257 0.6909 -58.87964 19.81450
## w1 - w3 18.29178 1.0000 -33.16718 69.75075
## w2 - w3 37.82435 0.1992 -11.75841 87.40712
print(out_16$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 462.2119 0.0000 *** 288.49560 635.9282
## w1 - w3 141.9628 0.3953 -85.22723 369.1528
## w2 - w3 -320.2491 0.0017 ** -539.15574 -101.3425
print(out_17$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 120.91225 1.0000 -250.4416 492.2661
## w1 - w3 204.77988 0.9235 -280.8849 690.4447
## w2 - w3 83.86763 1.0000 -384.0898 551.8250
print(out_18$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 -2.315477 1.0000 -39.54988 34.91892
## w1 - w3 25.776090 0.6032 -22.91988 74.47206
## w2 - w3 28.091567 0.4456 -18.82894 75.01208
print(out_19$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 529.7331 0.0000 *** 423.52655 635.9396
## w1 - w3 693.7035 0.0000 *** 554.80426 832.6027
## w2 - w3 163.9704 0.0107 * 30.13546 297.8054
print(out_20$comparison)
## difference pvalue signif. LCL UCL
## w1 - w2 920.3232 0.2248 -323.7992 2164.446
## w1 - w3 604.0234 1.0000 -1023.0674 2231.114
## w2 - w3 -316.2999 1.0000 -1884.0667 1251.467
#############################################################################################
###### Chi-square test with Benjamini and Hochberg false discovery-rate corrections #########
#############################################################################################
df3 = banco %>% select("record_id", "HMOSecretor", "HMOSecretor_1m", "HMOSecretor_3m")
headtail(df3)
## record_id HMOSecretor HMOSecretor_1m HMOSecretor_3m
## 1 1 NA NA NA
## 2 2 NA 1 NA
## 3 3 NA NA NA
## 146 146 NA NA NA
## 147 147 NA NA NA
## 148 148 NA NA NA
df3 = melt(df3, id="record_id")
headtail(df3)
## record_id variable value
## 1 1 HMOSecretor NA
## 2 2 HMOSecretor NA
## 3 3 HMOSecretor NA
## 442 146 HMOSecretor_3m NA
## 443 147 HMOSecretor_3m NA
## 444 148 HMOSecretor_3m NA
tab = table(df3$variable, df3$value)
tab
##
## 0 1
## HMOSecretor 6 37
## HMOSecretor_1m 6 49
## HMOSecretor_3m 2 19
chi_sq = chisq.test(tab) # simulate.p.value = T
## Warning in chisq.test(tab): Chi-squared approximation may be incorrect
chi_sq
## Pearson's Chi-squared test with tab
## X-squared = 0.3388, df = 2, p-value = 0.8442
chisqPostHoc(chi_sq, control = "BH")
## Warning in stats::chisq.test(tmp): Chi-squared approximation may be
## incorrect
## Warning in stats::chisq.test(tmp): Chi-squared approximation may be
## incorrect
## Adjusted p-values used the BH method.
## comparison raw.p adj.p
## 1 HMOSecretor vs. HMOSecretor_1m 0.8841 1
## 2 HMOSecretor vs. HMOSecretor_3m 0.9198 1
## 3 HMOSecretor_1m vs. HMOSecretor_3m 1.0000 1
#chisqPostHoc(chi_sq)
#chisq.test(tab[1,])
#chisq.test(tab[2,])
#chisq.test(tab[3,])
#################################################################
########### correlogram with the significance test ##############
#################################################################
bd_w1 = banco[,c(8,29,104,164:173,235:259)]
str(bd_w1)
## 'data.frame': 148 obs. of 38 variables:
## $ mother_age : num 36 29.4 25.8 27.7 23.1 ...
## $ parity : int NA NA NA NA NA 2 1 1 2 NA ...
## $ schooling : int 10 15 12 11 7 9 NA 11 5 12 ...
## $ total_energy : num 10287 5219 NA 3541 2811 ...
## $ total_protein : num 507 259 0 146 107 ...
## $ total_fat : num 434 203 0 129 102 ...
## $ total_carbohydrates: num 1120 603 0 460 374 ...
## $ total_fiber : num 43.1 42.5 0 47.5 22.3 ...
## $ saturated_fat : num 251.3 90.1 0 50.4 49.3 ...
## $ unsaturated_fat : num 131.4 75 0 55.5 38.4 ...
## $ calcium : num 14545 4258 0 1845 2379 ...
## $ iron : num 20.7 18.9 0 17.5 7.4 ...
## $ tryptophan : num 1.845 2.042 0 1.275 0.703 ...
## $ HMOSecretor : int NA NA NA NA NA NA NA NA NA NA ...
## $ Diversity_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ Evenness_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ 2'FL_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ 3FL_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNnT_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ 3'SL_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ DFLac_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ 6'SL_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNT_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNFPI_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNFPII_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNFPIII_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LSTb_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LSTc_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ DFLNT_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNH_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ DSLNT_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ FLNH_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ DFLNH_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ FDSLNH_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ DSLNH_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ SUM_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ Sia_7d : num NA NA NA NA NA NA NA NA NA NA ...
## $ Fuc_7d : num NA NA NA NA NA NA NA NA NA NA ...
bd_w2 = banco[,c(8,29,104,164:173,580:604)]
str(bd_w2)
## 'data.frame': 148 obs. of 38 variables:
## $ mother_age : num 36 29.4 25.8 27.7 23.1 ...
## $ parity : int NA NA NA NA NA 2 1 1 2 NA ...
## $ schooling : int 10 15 12 11 7 9 NA 11 5 12 ...
## $ total_energy : num 10287 5219 NA 3541 2811 ...
## $ total_protein : num 507 259 0 146 107 ...
## $ total_fat : num 434 203 0 129 102 ...
## $ total_carbohydrates: num 1120 603 0 460 374 ...
## $ total_fiber : num 43.1 42.5 0 47.5 22.3 ...
## $ saturated_fat : num 251.3 90.1 0 50.4 49.3 ...
## $ unsaturated_fat : num 131.4 75 0 55.5 38.4 ...
## $ calcium : num 14545 4258 0 1845 2379 ...
## $ iron : num 20.7 18.9 0 17.5 7.4 ...
## $ tryptophan : num 1.845 2.042 0 1.275 0.703 ...
## $ HMOSecretor_1m : int NA 1 NA NA NA 1 NA NA 1 NA ...
## $ Diversity_1m : num NA 3.42 NA NA NA ...
## $ Evenness_1m : num NA 0.18 NA NA NA ...
## $ 2'FL_1m : num NA 8227 NA NA NA ...
## $ 3FL_1m : num NA 290 NA NA NA ...
## $ LNnT_1m : num NA 172 NA NA NA ...
## $ 3'SL_1m : num NA 596 NA NA NA ...
## $ DFLac_1m : num NA 373 NA NA NA ...
## $ 6'SL_1m : num NA 374 NA NA NA ...
## $ LNT_1m : num NA 602 NA NA NA ...
## $ LNFPI_1m : num NA 2578 NA NA NA ...
## $ LNFPII_1m : num NA 608 NA NA NA ...
## $ LNFPIII_1m : num NA 27.1 NA NA NA ...
## $ LSTb_1m : num NA 66.1 NA NA NA ...
## $ LSTc_1m : num NA 235 NA NA NA ...
## $ DFLNT_1m : num NA 1633 NA NA NA ...
## $ LNH_1m : num NA 34.1 NA NA NA ...
## $ DSLNT_1m : num NA 83.8 NA NA NA ...
## $ FLNH_1m : num NA 56.7 NA NA NA ...
## $ DFLNH_1m : num NA 201 NA NA NA ...
## $ FDSLNH_1m : num NA 84.7 NA NA NA ...
## $ DSLNH_1m : num NA 161 NA NA NA ...
## $ SUM_1m : num NA 16401 NA NA NA ...
## $ Sia_1m : num NA 1929 NA NA NA ...
## $ Fuc_1m : num NA 16284 NA NA NA ...
bd_w3 = banco[,c(8,29,104,164:173,925:949)]
str(bd_w3)
## 'data.frame': 148 obs. of 38 variables:
## $ mother_age : num 36 29.4 25.8 27.7 23.1 ...
## $ parity : int NA NA NA NA NA 2 1 1 2 NA ...
## $ schooling : int 10 15 12 11 7 9 NA 11 5 12 ...
## $ total_energy : num 10287 5219 NA 3541 2811 ...
## $ total_protein : num 507 259 0 146 107 ...
## $ total_fat : num 434 203 0 129 102 ...
## $ total_carbohydrates: num 1120 603 0 460 374 ...
## $ total_fiber : num 43.1 42.5 0 47.5 22.3 ...
## $ saturated_fat : num 251.3 90.1 0 50.4 49.3 ...
## $ unsaturated_fat : num 131.4 75 0 55.5 38.4 ...
## $ calcium : num 14545 4258 0 1845 2379 ...
## $ iron : num 20.7 18.9 0 17.5 7.4 ...
## $ tryptophan : num 1.845 2.042 0 1.275 0.703 ...
## $ HMOSecretor_3m : int NA NA NA NA NA NA NA NA NA NA ...
## $ Diversity_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ Evenness_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ 2'FL_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ 3FL_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNnT_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ 3'SL_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ DFLac_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ 6'SL_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNT_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNFPI_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNFPII_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNFPIII_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LSTb_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LSTc_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ DFLNT_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ LNH_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ DSLNT_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ FLNH_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ DFLNH_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ FDSLNH_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ DSLNH_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ SUM_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ Sia_3m : num NA NA NA NA NA NA NA NA NA NA ...
## $ Fuc_3m : num NA NA NA NA NA NA NA NA NA NA ...
# colnames(bd_w1) = c()
res1 = cor.mtest(bd_w1, conf.level = .99, method="spearman", exact=F, use = "pairwise.complete.obs")
res2 = cor.mtest(bd_w2, conf.level = .99, method="spearman", exact=F, use = "pairwise.complete.obs")
res3 = cor.mtest(bd_w3, conf.level = .99, method="spearman", exact=F, use = "pairwise.complete.obs")
corrplot(cor(bd_w1,
method = "spearman", use = "pairwise.complete.obs"),
method = "circle", type = "upper", diag = F, tl.col = "black", number.cex = .3,
p.mat = res1$p, sig.level = 0.01)
corrplot(cor(bd_w2,
method = "spearman", use = "pairwise.complete.obs"),
method = "circle", type = "upper", diag = F, tl.col = "black", number.cex = .3,
p.mat = res2$p, sig.level = 0.01)
corrplot(cor(bd_w3,
method = "spearman", use = "pairwise.complete.obs"),
method = "circle", type = "upper", diag = F, tl.col = "black", number.cex = .3,
p.mat = res3$p, sig.level = 0.01)
## >> Exemplo de outra apresentação do corrplot
corrplot(cor(bd_w1,
method = "spearman", use = "pairwise.complete.obs"),
method = "circle", type = "upper", diag = F, tl.col = "black", number.cex = .3,
p.mat = res1$p, sig.level = 0.01,
insig = "blank") # leave blank on no significant coefficient
corrplot(cor(bd_w2,
method = "spearman", use = "pairwise.complete.obs"),
method = "circle", type = "upper", diag = F, tl.col = "black", number.cex = .3,
p.mat = res1$p, sig.level = 0.01,
insig = "blank") # leave blank on no significant coefficient
corrplot(cor(bd_w3,
method = "spearman", use = "pairwise.complete.obs"),
method = "circle", type = "upper", diag = F, tl.col = "black", number.cex = .3,
p.mat = res1$p, sig.level = 0.01,
insig = "blank") # leave blank on no significant coefficient
########################################
############## GGPLOTS #################
########################################
tb3 = banco %>% select("record_id", 238:256, 583:601, 928:946) # dim = 148 x 58 (=19 x 3)
tb3 = melt(tb3, id="record_id")
headtail(tb3)
## record_id variable value
## 1 1 2'FL_7d NA
## 2 2 2'FL_7d NA
## 3 3 2'FL_7d NA
## 8434 146 DSLNH_3m NA
## 8435 147 DSLNH_3m NA
## 8436 148 DSLNH_3m NA
tb3$onda=rep(c("w1","w2","w3"), each=2812) # 19 vars x 148 ind
tb3$onda=as.factor(tb3$onda)
headtail(tb3)
## record_id variable value onda
## 1 1 2'FL_7d NA w1
## 2 2 2'FL_7d NA w1
## 3 3 2'FL_7d NA w1
## 8434 146 DSLNH_3m NA w3
## 8435 147 DSLNH_3m NA w3
## 8436 148 DSLNH_3m NA w3
tb3$variable = as.factor(sub("[_].*", "", tb3$variable)) # levels(tb3$variable)
headtail(tb3)
## record_id variable value onda
## 1 1 2'FL NA w1
## 2 2 2'FL NA w1
## 3 3 2'FL NA w1
## 8434 146 DSLNH NA w3
## 8435 147 DSLNH NA w3
## 8436 148 DSLNH NA w3
sum_hmo = Rmisc::summarySE(tb3, measurevar="value", groupvars=c("onda","variable"),
na.rm = T)
headtail(sum_hmo)
## onda variable N value sd se ci
## 1 w1 2'FL 43 4575.30384 2394.26944 365.122666 736.84737
## 2 w1 3'SL 43 342.45873 187.23204 28.552619 57.62152
## 3 w1 3FL 43 235.13986 116.77472 17.807977 35.93795
## 55 w3 LNT 21 1412.98045 674.98590 147.293998 307.24990
## 56 w3 LSTb 21 94.99957 38.70427 8.445964 17.61797
## 57 w3 LSTc 21 131.57653 78.46418 17.122288 35.71647
# ABSOLUTE GGPLOT
library(RColorBrewer)
colourCount = length(unique(sum_hmo$variable))
getPalette = colorRampPalette(brewer.pal(12, "Paired")) # "Set3" "Accent" "Set1"
#dev.new()
ggplot(data=sum_hmo, aes(x=onda, y=value, group=variable, fill=variable)) +
geom_bar(stat="identity", position=position_stack(reverse=T)) +
scale_fill_manual(values=getPalette(colourCount),
guide = guide_legend(reverse=T)) +
labs(title="HMO profile by postpartum period",
x="Postpartum period",
y="HMO concentration (nmol/mL)") +
scale_x_discrete(labels=c("w1"="1-7 days","w2"="30-45 days","w3"="3 months")) +
theme_minimal() +
theme(panel.grid = element_blank(),
legend.position="right",
legend.title=element_blank(),
text = element_text(size=13),
axis.text=element_text(size = 13))
# RELATIVE GGPLOT
#dev.new()
ggplot(data=sum_hmo, aes(x=onda, y=value, group=variable, fill=variable)) +
geom_bar(stat="identity", position=position_fill(reverse = TRUE)) +
scale_y_continuous(labels=scales::percent) +
scale_fill_manual(values=getPalette(colourCount),
guide = guide_legend(reverse=T)) +
labs(title="HMO profile by postpartum period",
x="Postpartum period",
y="Relative HMO concentration (% of total)") +
scale_x_discrete(labels=c("w1"="1-7 days","w2"="30-45 days","w3"="3 months")) +
theme_minimal() +
theme(panel.grid = element_blank(),
legend.position="right",
legend.title=element_blank(),
text = element_text(size=13),
axis.text=element_text(size = 13))
#############
#### FIM ####
#############
#require(graphics) # library(help = "graphics")
#dev.new()
#pairs(bd_w1[17:27]) # matrix of scatterplots
#pairs(bd_w1[28:36]) # matrix of scatterplots
#dev.off()
#rm(list=ls())
#q()