Input data and packages
library(readxl)
## Warning: package 'readxl' was built under R version 4.0.3
Data <- read_excel("C:/Users/Admin/Desktop/R/Data.xlsx",
col_types = c("text", "text", "text",
"numeric", "numeric", "numeric",
"numeric", "numeric","numeric"))
require(multcomp)
## Loading required package: multcomp
## Warning: package 'multcomp' was built under R version 4.0.3
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.0.3
## Loading required package: survival
## Loading required package: TH.data
## Warning: package 'TH.data' was built under R version 4.0.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.0.3
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
require(lsmeans)
## Loading required package: lsmeans
## Warning: package 'lsmeans' was built under R version 4.0.3
## Loading required package: emmeans
## Warning: package 'emmeans' was built under R version 4.0.3
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
require(ggpubr)
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.0.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.3
require(ggplot2)
attach(Data)
Define factors
Data$`Feeding level`=as.factor(Data$`Feeding level`)
Data$Treatment=as.factor(Data$Treatment)
str(Data)
## tibble [24 x 9] (S3: tbl_df/tbl/data.frame)
## $ Feeding level : Factor w/ 2 levels "Less","More": 1 1 1 2 2 2 1 1 1 2 ...
## $ Treatment : Factor w/ 4 levels "Base","High carb",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ Gender : chr [1:24] "Male" "Male" "Male" "Male" ...
## $ Body Weight (g) : num [1:24] 2506 1433 1555 2806 2530 ...
## $ Liver (g) : num [1:24] 30.6 27.6 22.6 56.4 44.8 ...
## $ Viscera (g) : num [1:24] 68.9 72 69.7 147.2 102.5 ...
## $ Fillet (g) : num [1:24] 442 492 440 987 796 ...
## $ Abdomial fat (g): num [1:24] 3.84 6.27 6.62 26.59 15.62 ...
## $ Weight gain (g) : num [1:24] 142 121 84.2 74.6 111 131 129 177 163 191 ...
Visualiation data, maybe no interaction between Treatment and Feeding level, the more feeding level led to the more of BW
plot=ggline(Data,x="Treatment",y="Body Weight (g)",color="Feeding level",add=c("mean_se","dotplot"),palette = c("black","pink"))
plot
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
RUN ANOVA, checking interaction
model=lm(`Body Weight (g)`~Treatment*`Feeding level`,data=Data)
anova(model)
## Analysis of Variance Table
##
## Response: Body Weight (g)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 1265961 421987 6.8517 0.003513 **
## `Feeding level` 1 4984082 4984082 80.9249 1.172e-07 ***
## Treatment:`Feeding level` 3 108185 36062 0.5855 0.633164
## Residuals 16 985424 61589
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No interaction (p=0.633), then checking in each factor, the result shows that each factor significantly affected to Body weight
model_no_inter=lm(`Body Weight (g)`~Treatment+`Feeding level`,data=Data)
anova(model_no_inter)
## Analysis of Variance Table
##
## Response: Body Weight (g)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 1265961 421987 7.3315 0.001849 **
## `Feeding level` 1 4984082 4984082 86.5918 1.656e-08 ***
## Residuals 19 1093609 57558
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Assumption 1: All samples are indepentdent (no repeated measured in 1 individual)
Assumption 2: Denpendent variable is continuous (all denpendent variables are already continous)
Assumption 3: Normal distribution (shapiro.test)
Assumption 4: Homogeneity of variances (bartlett.test)
Data visualization (No interaction => boxplot)
d=ggplot(Data,aes(x=Treatment,y=`Body Weight (g)`,fill=`Feeding level`))
d+geom_boxplot()+theme_bw()+theme_classic()
Posthoc
posthoc=lsmeans(model_no_inter,pairwise~Treatment,adjust="Tukey")
cld(posthoc[[1]],alpha=0.05,Letters=letters)
## Treatment lsmean SE df lower.CL upper.CL .group
## High carb 1644 97.9 19 1439 1849 a
## Mix 1699 97.9 19 1494 1904 a
## High fat 2003 97.9 19 1798 2208 ab
## Base 2208 97.9 19 2003 2413 b
##
## Results are averaged over the levels of: Feeding level
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc1=lsmeans(model_no_inter,pairwise~`Feeding level`,adjust="Tukey")
cld(posthoc1[[1]],alpha=0.05,Letters=letters)
## Feeding level lsmean SE df lower.CL upper.CL .group
## Less 1433 69.3 19 1288 1578 a
## More 2344 69.3 19 2199 2489 b
##
## Results are averaged over the levels of: Treatment
## Confidence level used: 0.95
## significance level used: alpha = 0.05