Data summary
#change column names
oldnames = c("Initial_ind_weight","Final_ind_weight")
newnames = c("Initial_weight","Final_individual_weight")
colnames(data)[colnames(data) %in% oldnames] <- newnames
attach(data)
#summary of data
library(psych)
str(data)## Classes 'tbl_df', 'tbl' and 'data.frame': 144 obs. of 22 variables:
## $ Tank : num 69 141 128 115 110 38 137 140 62 143 ...
## $ Label : chr "F" "F" "X" "F" ...
## $ CP : num 37 37 42 37 36 40 42 40 42 47 ...
## $ Lysine : num 3 3 3.3 3 2.7 3 3.3 3 3.3 3.6 ...
## $ M_C : num 1.71 1.71 1.9 1.71 1.54 1.71 1.9 1.71 1.9 2.05 ...
## $ Treatment : chr "3.0% lysine. -6% CP" "3.0% lysine. -6% CP" "3.3% lysine. -6% CP" "3.0% lysine. -6% CP" ...
## $ Initial_count : num 7 7 7 7 7 7 7 7 7 7 ...
## $ Initial_group_weight : num 36.9 33.5 37.7 34.8 36.6 34.4 35.7 36.3 36.8 34.8 ...
## $ Initial_weight : num 5.27 4.79 5.39 4.97 5.23 ...
## $ Final_count : num 7 7 7 7 7 7 7 7 7 6 ...
## $ Final_group_weight : num 158 154 146 144 144 ...
## $ Final_individual_weight: num 22.6 22 20.9 20.6 20.6 ...
## $ MIA : num NA NA NA NA NA NA NA NA NA NA ...
## $ Feed_orig : num 187 187 187 187 187 ...
## $ Feed : num 26.7 26.7 26.7 26.7 26.7 ...
## $ Days : num 43 43 43 43 43 43 43 43 43 43 ...
## $ Weight_gain : num 328 361 287 314 293 ...
## $ Weight_gain_orig : num 17.3 17.3 15.5 15.6 15.3 ...
## $ ADG : num 0.525 0.513 0.485 0.479 0.478 ...
## $ Biomass_gain : num 121 121 108 109 107 ...
## $ FCR : num 1.54 1.55 1.73 1.71 1.74 ...
## $ Survival : num 100 100 100 100 100 ...#Mean, Sum of each Treatment
data <- data[,c(2,16,18,19)]
#Subset data
data$Label[data$Label == "Ctrl 1"] <- "Ctrl1"
data$Label[data$Label == "Ctrl 2"] <- "Ctrl2"
data$Label <- as.factor(data$Label)
data %>% group_by(Label) %>% summarise_each(funs(mean, sd))## # A tibble: 18 x 7
## Label Days_mean Weight_gain_orig_m~ ADG_mean Days_sd Weight_gain_orig~ ADG_sd
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A 42.5 13.2 0.430 0.535 1.54 0.0356
## 2 B 42.5 13.8 0.442 0.535 0.897 0.0224
## 3 C 42.5 13.4 0.434 0.535 1.43 0.0370
## 4 Ctrl1 42.5 11.0 0.379 0.535 1.40 0.0333
## 5 Ctrl2 42.5 10.1 0.359 0.535 0.838 0.0198
## 6 D 42.5 13.1 0.428 0.535 0.955 0.0230
## 7 E 42.5 14.1 0.451 0.535 1.09 0.0249
## 8 F 42.5 15.3 0.481 0.535 1.35 0.0271
## 9 G 42.5 12.8 0.417 0.535 1.56 0.0373
## 10 H 42.5 12.6 0.417 0.535 0.891 0.0230
## 11 K 42.5 10.4 0.364 0.535 3.81 0.0881
## 12 L 42.5 13.3 0.434 0.535 1.76 0.0444
## 13 P 42.5 12.7 0.419 0.535 1.37 0.0315
## 14 Q 42.5 7.18 0.287 0.535 1.03 0.0302
## 15 R 42.5 11.9 0.402 0.535 1.01 0.0259
## 16 X 42.5 14.4 0.459 0.535 1.27 0.0283
## 17 Y 42.5 11.9 0.399 0.535 1.93 0.0490
## 18 Z 42.5 3.93 0.213 0.535 1.28 0.0298# plotmeans
plotmeans(ADG ~ Label, data = data, frame = FALSE,
xlab = "Treatment", ylab = "ADG",
main="Mean Plot with 95% CI") 
Assumption Test
Homogeneity of variances Test
# 1. Homogeneity of variances
#leveneTest
library(car)
levene <- vector(mode = "list", length = 3)
anova <- vector(mode = "list", length = 3)
for (i in 2:4){
levene[[i-1]] <-leveneTest(data[[i]]~data[[1]],data=data, center=median) #default center is mean
cat(paste("* ", colnames(data)[i]),'\n')
print(levene[[i-1]])
anova[[i-1]] <-aov(data[[i]]~data[[1]],data=data)
plot(anova[[i-1]], 1)
}## * Days
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 1 0.4361
## 49
## * Weight_gain_orig
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 1.9361 0.09355 .
## 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## * ADG
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 1.8781 0.1035
## 49
H0 is that all variances are equal. The test reveals a p-value greater than 0.05, indicating that there is no significant difference between the group variances in location.
Normality Test
# 2. Normality
# Run Shapiro-Wilk test
normtest <- vector(mode = "list", length = 3)
anova <- vector(mode = "list", length = 3)
for (i in 2:4){
normtest[[i-1]] <-shapiro.test(data[[i]])
cat(paste("* ", colnames(data)[i]),'\n')
print(normtest[[i-1]]$p.value)
anova[[i-1]] <-aov(data[[i]]~data[[1]],data=data)
plot(anova[[i-1]], 2)
}## * Days
## [1] 1.628755e-10
## * Weight_gain_orig
## [1] 0.000283818
## * ADG
## [1] 0.0004340944
Data transformation
#square root
st <- function(p) { sqrt(p) }
stnorm <- vector(mode = "list", length = 5)
for (i in 2:4){
stnorm[[i-1]] <-shapiro.test(st(data[[i]]))
cat(paste("* ", colnames(data)[i]),'\n')
print(stnorm[[i-1]]$p.value)
}## * Days
## [1] 1.628755e-10
## * Weight_gain_orig
## [1] 8.135325e-07
## * ADG
## [1] 1.548732e-05#1/square root
re_sqrt <- function(p) { 1/sqrt(p) }
#logit
logitTransform <- function(p) { log(p/(1-p)) }
logitnorm <- vector(mode = "list", length = 5)
for (i in 2:4){
logitnorm[[i-1]] <-shapiro.test(st(data[[i]]))
cat(paste("* ", colnames(data)[i]),'\n')
print(logitnorm[[i-1]]$p.value)
}## * Days
## [1] 1.628755e-10
## * Weight_gain_orig
## [1] 8.135325e-07
## * ADG
## [1] 1.548732e-05#log
logTransform <- function(p) { log(p) }
lognorm <- vector(mode = "list", length = 5)
for (i in 2:4){
lognorm[[i-1]] <-shapiro.test(st(data[[i]]))
cat(paste("* ", colnames(data)[i]),'\n')
print(lognorm[[i-1]]$p.value)
}## * Days
## [1] 1.628755e-10
## * Weight_gain_orig
## [1] 8.135325e-07
## * ADG
## [1] 1.548732e-05parametric ANOVA
# ANOVA Test - Parametric
anova <- vector(mode = "list", length = 2)
for (i in 2:4){
anova[[i-1]] <-aov(data[[i]]~data[[1]],data=data)
cat(paste("* ", colnames(data)[i]),'\n')
print(summary(anova[[i-1]]))
}## * Days
## Df Sum Sq Mean Sq F value Pr(>F)
## data[[1]] 6 0 0.0000 0 1
## Residuals 49 14 0.2857
## * Weight_gain_orig
## Df Sum Sq Mean Sq F value Pr(>F)
## data[[1]] 6 96.42 16.070 4.839 0.000597 ***
## Residuals 49 162.72 3.321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## * ADG
## Df Sum Sq Mean Sq F value Pr(>F)
## data[[1]] 6 0.05184 0.008640 4.767 0.000673 ***
## Residuals 49 0.08882 0.001813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Post-hoc Analysis: letter grouping (1)
#Tukey HSD (Tukey Honest Significant Differences)
library(agricolae)
anova <- vector(mode = "list", length = 2)
for (i in 2:4){
anova[[i-1]] <-aov(data[[i]]~data[[1]],data=data)
cat(paste("* ", colnames(data)[i]),'\n')
print(TukeyHSD(anova[[i-1]]))
print(HSD.test(anova[[i-1]],"data[[1]]", group=TRUE)$'groups')
}## * Days
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = data[[i]] ~ data[[1]], data = data)
##
## $`data[[1]]`
## diff lwr upr p adj
## B-A 7.105427e-15 -0.8215703 0.8215703 1
## Ctrl1-A 0.000000e+00 -0.8215703 0.8215703 1
## Ctrl2-A 0.000000e+00 -0.8215703 0.8215703 1
## K-A 0.000000e+00 -0.8215703 0.8215703 1
## P-A 0.000000e+00 -0.8215703 0.8215703 1
## R-A 0.000000e+00 -0.8215703 0.8215703 1
## Ctrl1-B -7.105427e-15 -0.8215703 0.8215703 1
## Ctrl2-B -7.105427e-15 -0.8215703 0.8215703 1
## K-B -7.105427e-15 -0.8215703 0.8215703 1
## P-B -7.105427e-15 -0.8215703 0.8215703 1
## R-B -7.105427e-15 -0.8215703 0.8215703 1
## Ctrl2-Ctrl1 0.000000e+00 -0.8215703 0.8215703 1
## K-Ctrl1 0.000000e+00 -0.8215703 0.8215703 1
## P-Ctrl1 0.000000e+00 -0.8215703 0.8215703 1
## R-Ctrl1 0.000000e+00 -0.8215703 0.8215703 1
## K-Ctrl2 0.000000e+00 -0.8215703 0.8215703 1
## P-Ctrl2 0.000000e+00 -0.8215703 0.8215703 1
## R-Ctrl2 0.000000e+00 -0.8215703 0.8215703 1
## P-K 0.000000e+00 -0.8215703 0.8215703 1
## R-K 0.000000e+00 -0.8215703 0.8215703 1
## R-P 0.000000e+00 -0.8215703 0.8215703 1
##
## data[[i]] groups
## A 42.5 a
## B 42.5 a
## Ctrl1 42.5 a
## Ctrl2 42.5 a
## K 42.5 a
## P 42.5 a
## R 42.5 a
## * Weight_gain_orig
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = data[[i]] ~ data[[1]], data = data)
##
## $`data[[1]]`
## diff lwr upr p adj
## B-A 0.5922619 -2.2086373 3.39316113 0.9946228
## Ctrl1-A -2.1875000 -4.9883992 0.61339922 0.2202393
## Ctrl2-A -3.0523810 -5.8532802 -0.25148173 0.0244551
## K-A -2.7795238 -5.5804230 0.02137542 0.0529979
## P-A -0.4437500 -3.2446492 2.35714922 0.9989175
## R-A -1.2592262 -4.0601254 1.54167303 0.8085569
## Ctrl1-B -2.7797619 -5.5806611 0.02113732 0.0529636
## Ctrl2-B -3.6446429 -6.4455421 -0.84374363 0.0037718
## K-B -3.3717857 -6.1726849 -0.57088649 0.0091825
## P-B -1.0360119 -3.8369111 1.76488732 0.9133274
## R-B -1.8514881 -4.6523873 0.94941113 0.4084777
## Ctrl2-Ctrl1 -0.8648810 -3.6657802 1.93601827 0.9622293
## K-Ctrl1 -0.5920238 -3.3929230 2.20887542 0.9946345
## P-Ctrl1 1.7437500 -1.0571492 4.54464922 0.4813162
## R-Ctrl1 0.9282738 -1.8726254 3.72917303 0.9471997
## K-Ctrl2 0.2728571 -2.5280421 3.07375637 0.9999350
## P-Ctrl2 2.6086310 -0.1922683 5.40953018 0.0831020
## R-Ctrl2 1.7931548 -1.0077445 4.59405399 0.4473859
## P-K 2.3357738 -0.4651254 5.13667303 0.1599783
## R-K 1.5202976 -1.2806016 4.32119684 0.6395132
## R-P -0.8154762 -3.6163754 1.98542303 0.9716475
##
## data[[i]] groups
## B 13.76964 a
## A 13.17738 ab
## P 12.73363 abc
## R 11.91815 abc
## Ctrl1 10.98988 abc
## K 10.39786 bc
## Ctrl2 10.12500 c
## * ADG
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = data[[i]] ~ data[[1]], data = data)
##
## $`data[[1]]`
## diff lwr upr p adj
## B-A 0.012034290 -0.053404023 0.0774726044 0.9975006
## Ctrl1-A -0.051217667 -0.116655981 0.0142206466 0.2180990
## Ctrl2-A -0.071376826 -0.136815140 -0.0059385120 0.0242612
## K-A -0.066024561 -0.131462875 -0.0005862471 0.0466715
## P-A -0.011144564 -0.076582878 0.0542937495 0.9983725
## R-A -0.028319128 -0.093757442 0.0371191859 0.8343722
## Ctrl1-B -0.063251958 -0.128690272 0.0021863561 0.0643140
## Ctrl2-B -0.083411116 -0.148849430 -0.0179728025 0.0048294
## K-B -0.078058851 -0.143497165 -0.0126205375 0.0101203
## P-B -0.023178855 -0.088617169 0.0422594590 0.9284445
## R-B -0.040353418 -0.105791732 0.0250848954 0.4928354
## Ctrl2-Ctrl1 -0.020159159 -0.085597473 0.0452791553 0.9626536
## K-Ctrl1 -0.014806894 -0.080245208 0.0506314202 0.9922603
## P-Ctrl1 0.040073103 -0.025365211 0.1055114168 0.5012350
## R-Ctrl1 0.022898539 -0.042539775 0.0883368532 0.9322578
## K-Ctrl2 0.005352265 -0.060086049 0.0707905788 0.9999768
## P-Ctrl2 0.060232262 -0.005206052 0.1256705754 0.0897870
## R-Ctrl2 0.043057698 -0.022380616 0.1084960118 0.4140760
## P-K 0.054879997 -0.010558317 0.1203183105 0.1552863
## R-K 0.037705433 -0.027732881 0.1031437468 0.5730969
## R-P -0.017174564 -0.082612878 0.0482637503 0.9831781
##
## data[[i]] groups
## B 0.4422926 a
## A 0.4302583 a
## P 0.4191137 ab
## R 0.4019391 ab
## Ctrl1 0.3790406 ab
## K 0.3642337 b
## Ctrl2 0.3588814 bPost-hoc Analysis: letter grouping (2)
library(multcomp)
Tukey <- vector(mode = "list", length = 3)
for (i in 2:4){
anova[[i-1]] <-aov(data[[i]]~Label,data=data)
cat(paste("* ", colnames(data)[i]),'\n')
print(cld(summary(glht(anova[[i-1]], linfct = mcp(Label = "Tukey")))))
old.par <- par(mai=c(1,1,1.25,1), no.readonly = TRUE)
plot(cld(summary(glht(anova[[i-1]], linfct = mcp(Label = "Tukey")))))
par(old.par)
}## * Days
## A B Ctrl1 Ctrl2 K P R
## "a" "a" "a" "a" "a" "a" "a"
## * Weight_gain_orig
## A B Ctrl1 Ctrl2 K P R
## "bc" "c" "ac" "a" "ab" "ac" "ac"
## * ADG
## A B Ctrl1 Ctrl2 K P R
## "b" "b" "ab" "a" "a" "ab" "ab"
non-parametric ANOVA
Kruskal-wallis test is an analysis of variancc performed on ranks. If the distributons among group are similars, Kruskal-Wallis test can test for a difference in medians. When the distributions are significantly different, rejection of null hypothesis can be happened by lack of independence. If we found differences among different levels of a group, post-hoc analysis can be performed to determine which levels of the independent variable differ from each other level.
kruskal <- vector(mode = "list", length = 2)
for (i in 2:4){
kruskal[[i-1]] <-kruskal.test(data[[i]]~data[[1]],data=data)
cat(paste("* ", colnames(data)[i]),'\n')
print(kruskal[[i-1]])
}## * Days
##
## Kruskal-Wallis rank sum test
##
## data: data[[i]] by data[[1]]
## Kruskal-Wallis chi-squared = 0, df = 6, p-value = 1
##
## * Weight_gain_orig
##
## Kruskal-Wallis rank sum test
##
## data: data[[i]] by data[[1]]
## Kruskal-Wallis chi-squared = 27.844, df = 6, p-value = 0.0001005
##
## * ADG
##
## Kruskal-Wallis rank sum test
##
## data: data[[i]] by data[[1]]
## Kruskal-Wallis chi-squared = 26.3, df = 6, p-value = 0.0001957post-hoc tests(Dunn’s test)
Dunn’s test can even apply to unequal number of observations in each level of the group.
To be free from type I error, p-value adjustment needs to be made. The options can be; “sidak”, “holm”, “hs”, “hochberg”, “bh”, “by”, “bonferroni”
library(rcompanion)
library(FSA)
Dunn <- vector(mode = "list", length = 2)
for (i in c(3)){
Dunn[[i-1]] <-dunnTest(data[[i]]~data[[1]],data=data,method="bh")
PT = Dunn[[i-1]]$res
print(PT)
print(cldList(comparison = Dunn[[i-1]]$res$Comparison,
p.value = Dunn[[i-1]]$res$P.adj,
threshold = 0.05))
}## Comparison Z P.unadj P.adj
## 1 A - B -0.7204387 4.712549e-01 0.5497974006
## 2 A - Ctrl1 2.4985428 1.247051e-02 0.0436467779
## 3 B - Ctrl1 3.2189815 1.286468e-03 0.0090052739
## 4 A - Ctrl2 3.6788360 2.343008e-04 0.0024601587
## 5 B - Ctrl2 4.3992748 1.086132e-05 0.0002280878
## 6 Ctrl1 - Ctrl2 1.1802932 2.378836e-01 0.3330370546
## 7 A - K 1.9007320 5.733713e-02 0.1204079753
## 8 B - K 2.6211707 8.762838e-03 0.0368039193
## 9 Ctrl1 - K -0.5978109 5.499661e-01 0.6078573195
## 10 Ctrl2 - K -1.7781041 7.538676e-02 0.1439201815
## 11 A - P 0.5058400 6.129690e-01 0.6436174575
## 12 B - P 1.2262787 2.200938e-01 0.3301407329
## 13 Ctrl1 - P -1.9927029 4.629399e-02 0.1080193211
## 14 Ctrl2 - P -3.1729961 1.508746e-03 0.0079209144
## 15 K - P -1.3948920 1.630484e-01 0.2633858943
## 16 A - R 1.5788338 1.143742e-01 0.2001548217
## 17 B - R 2.2992725 2.148947e-02 0.0644684072
## 18 Ctrl1 - R -0.9197090 3.577248e-01 0.4418953948
## 19 Ctrl2 - R -2.1000022 3.572864e-02 0.0937876912
## 20 K - R -0.3218982 7.475299e-01 0.7475298523
## 21 P - R 1.0729938 2.832739e-01 0.3717969651
## Group Letter MonoLetter
## 1 A ab ab
## 2 B a a
## 3 Ctrl1 cd cd
## 4 Ctrl2 c c
## 5 K bcd bcd
## 6 P abd ab d
## 7 R abcd abcdpost-hoc tests(Nemenyi test)
Nemenyi test is one of the post-hoc method of Kruskal-walli. It is not appropriate for groups with unequal numbers of observations. “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”
#Pairwise comparisons using Tukey and Kramer (Nemenyi) test with Tukey-Dist approximation for independent samples
library(PMCMRplus)
data1 <- data
nemenyi <- vector(mode = "list", length = 3)
for (i in 2:4){
cat(paste("* ", colnames(data)[i]),'\n')
nemenyi[[i-1]] <-kwAllPairsNemenyiTest(data1[[i]] ~ data1$Label,data = data1)
print(summary(nemenyi[[i-1]]))
plot(nemenyi[[i-1]])
}## * Days
## q value Pr(>|q|)
## B - A == 0 0 1
## Ctrl1 - A == 0 0 1
## Ctrl2 - A == 0 0 1
## K - A == 0 0 1
## P - A == 0 0 1
## R - A == 0 0 1
## Ctrl1 - B == 0 0 1
## Ctrl2 - B == 0 0 1
## K - B == 0 0 1
## P - B == 0 0 1
## R - B == 0 0 1
## Ctrl2 - Ctrl1 == 0 0 1
## K - Ctrl1 == 0 0 1
## P - Ctrl1 == 0 0 1
## R - Ctrl1 == 0 0 1
## K - Ctrl2 == 0 0 1
## P - Ctrl2 == 0 0 1
## R - Ctrl2 == 0 0 1
## P - K == 0 0 1
## R - K == 0 0 1
## R - P == 0 0 1
## A B Ctrl1 Ctrl2 K P
## B 1 - - - - -
## Ctrl1 1 1 - - - -
## Ctrl2 1 1 1 - - -
## K 1 1 1 1 - -
## P 1 1 1 1 1 -
## R 1 1 1 1 1 1
## * Weight_gain_orig
## q value Pr(>|q|)
## B - A == 0 1.019 0.99140400
## Ctrl1 - A == 0 3.533 0.15963880
## Ctrl2 - A == 0 5.203 0.00439149 **
## K - A == 0 2.688 0.47973316
## P - A == 0 0.715 0.99878562
## R - A == 0 2.233 0.69602080
## Ctrl1 - B == 0 4.552 0.02189595 *
## Ctrl2 - B == 0 6.222 0.00021946 ***
## K - B == 0 3.707 0.11953179
## P - B == 0 1.734 0.88415088
## R - B == 0 3.252 0.24416080
## Ctrl2 - Ctrl1 == 0 1.669 0.90177705
## K - Ctrl1 == 0 0.845 0.99689614
## P - Ctrl1 == 0 2.818 0.41922734
## R - Ctrl1 == 0 1.301 0.96955349
## K - Ctrl2 == 0 2.515 0.56295266
## P - Ctrl2 == 0 4.487 0.02535673 *
## R - Ctrl2 == 0 2.970 0.35251051
## P - K == 0 1.973 0.80485595
## R - K == 0 0.455 0.99991117
## R - P == 0 1.517 0.93608097
## A B Ctrl1 Ctrl2 K P
## B 0.99140 - - - - -
## Ctrl1 0.15964 0.02190 - - - -
## Ctrl2 0.00439 0.00022 0.90178 - - -
## K 0.47973 0.11953 0.99690 0.56295 - -
## P 0.99879 0.88415 0.41923 0.02536 0.80486 -
## R 0.69602 0.24416 0.96955 0.35251 0.99991 0.93608
## * ADG
## q value Pr(>|q|)
## B - A == 0 1.041 0.9903806
## Ctrl1 - A == 0 3.501 0.1681449
## Ctrl2 - A == 0 4.964 0.0081467 **
## K - A == 0 2.645 0.5003663
## P - A == 0 0.650 0.9992950
## R - A == 0 2.027 0.7838614
## Ctrl1 - B == 0 4.541 0.0224421 *
## Ctrl2 - B == 0 6.005 0.0004351 ***
## K - B == 0 3.685 0.1240752
## P - B == 0 1.691 0.8960978
## R - B == 0 3.067 0.3124934
## Ctrl2 - Ctrl1 == 0 1.463 0.9460811
## K - Ctrl1 == 0 0.856 0.9966689
## P - Ctrl1 == 0 2.851 0.4045286
## R - Ctrl1 == 0 1.474 0.9441725
## K - Ctrl2 == 0 2.320 0.6560971
## P - Ctrl2 == 0 4.314 0.0370225 *
## R - Ctrl2 == 0 2.937 0.3663807
## P - K == 0 1.994 0.7965789
## R - K == 0 0.618 0.9994749
## R - P == 0 1.377 0.9597638
## A B Ctrl1 Ctrl2 K P
## B 0.99038 - - - - -
## Ctrl1 0.16814 0.02244 - - - -
## Ctrl2 0.00815 0.00044 0.94608 - - -
## K 0.50037 0.12408 0.99667 0.65610 - -
## P 0.99929 0.89610 0.40453 0.03702 0.79658 -
## R 0.78386 0.31249 0.94417 0.36638 0.99947 0.95976