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",
sheet = "2 way interaction", col_types = c("text",
"text", "numeric"))
View(Data)
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(car)
## Loading required package: car
## Warning: package 'car' was built under R version 4.0.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.0.3
attach(Data)
Define factors
Data$Treatment=as.factor(Data$Treatment)
Data$Gender=as.factor(Data$Gender)
str(Data)
## tibble [24 x 3] (S3: tbl_df/tbl/data.frame)
## $ Gender : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
## $ Treatment : Factor w/ 2 levels "Base","High carb": 1 1 1 1 1 1 2 2 2 2 ...
## $ Weight gain (g): num [1:24] 142 121 84.2 74.6 111 131 129 177 163 191 ...
Visualiation data, maybe there have the interaction between Treatment and Gender.
plot=ggline(Data,x="Treatment",y="Weight gain (g)",color="Gender",add=c("mean_se","dotplot"),palette = c("black","pink"))
plot
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
RUN ANOVA, checking interaction => There have interaction (p=0.004)
model=lm(`Weight gain (g)`~Gender*Treatment,data=Data)
anova(model)
## Analysis of Variance Table
##
## Response: Weight gain (g)
## Df Sum Sq Mean Sq F value Pr(>F)
## Gender 1 27696 27696 30.467 2.107e-05 ***
## Treatment 1 57145 57145 62.861 1.339e-07 ***
## Gender:Treatment 1 9389 9389 10.328 0.004357 **
## Residuals 20 18181 909
## ---
## 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
anova_residual=residuals(object = model)
shapiro.test(anova_residual)
##
## Shapiro-Wilk normality test
##
## data: anova_residual
## W = 0.95521, p-value = 0.3498
plot(model,2) # mean:qqnorm(`Weight gain (g)`)+ qqline(`Weight gain (g)`)
Assumption 4: Homogeneity of variances
plot(model,1)
leveneTest(`Weight gain (g)`~Treatment*Gender,data=Data)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.1476 0.93
## 20
Posthoc
posthoc=lsmeans(model,pairwise~Treatment*Gender,adjust="Tukey")
cld(posthoc[[1]],alpha=0.05,Letters=letters)
## Treatment Gender lsmean SE df lower.CL upper.CL .group
## Base Male 111 12.3 20 85 136 a
## Base Female 139 12.3 20 113 165 ab
## High carb Male 169 12.3 20 143 194 b
## High carb Female 276 12.3 20 250 302 c
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05