Analysis of variance

#reading data

library(HH)
library(gplots)
library(multcomp)
library(MASS)
library(tidyverse)
library(readxl)
library(viridis)
library(hrbrthemes)
 fc <- "c:\\Users\\soroush\\Desktop\\kimia\\Proximates workbook.xlsx"
 dat <- read_excel(fc,sheet=1)
 data <- dat[-1,]
 data <- data[c(-13,-26,-39,-52,-65),]

#change strings to factors

str(data)

## Classes 'tbl_df', 'tbl' and 'data.frame':    60 obs. of  4 variables:
##  $ Time_point: chr  "T6" "T6" "T6" "T6" ...
##  $ % Fat     : num  4.1 3.48 3.04 2.3 3.49 ...
##  $ % Protein : num  20.6 20.4 21 20.6 20 ...
##  $ group     : chr  "FD-/RF-" "FD+/RF+" "FD-/RF+" "FD-/RF-" ...

data$Time_point <- factor(data$Time_point,levels = c("T6","T7","T8","T9","T10"))
levels(data$Time_point)[levels(data$Time_point)=="T6"] <- "initial"
levels(data$Time_point)[levels(data$Time_point)=="T7"] <- "1st fasting"
levels(data$Time_point)[levels(data$Time_point)=="T8"] <- "1st refeeding"
levels(data$Time_point)[levels(data$Time_point)=="T9"] <- "2nd fasting"
levels(data$Time_point)[levels(data$Time_point)=="T10"] <- "2nd refeeding"
data$group <- factor(data$group)
table(data$group)

## 
## FD-/RF- FD-/RF+ FD+/RF- FD+/RF+ 
##      15      15      15      15

table(data$Time_point)

## 
##       initial   1st fasting 1st refeeding   2nd fasting 2nd refeeding 
##            12            12            12            12            12

#two way manova model

attach(data)
names(data)[2] <- "fat"
names(data)[3] <- "protein"

y <-  cbind(data$fat,data$protein)

fit <- manova(y ~ Time_point*group , data=data)
summary(fit)

##                  Df  Pillai approx F num Df den Df    Pr(>F)    
## Time_point        4 0.91817   8.4873      8     80 2.754e-08 ***
## group             3 0.11403   0.8062      6     80    0.5681    
## Time_point:group 12 0.42369   0.8959     24     80    0.6062    
## Residuals        40                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(fit)

##  Response 1 :
##                  Df Sum Sq Mean Sq F value    Pr(>F)    
## Time_point        4 27.018  6.7545  7.3834 0.0001495 ***
## group             3  0.081  0.0271  0.0297 0.9930074    
## Time_point:group 12  4.260  0.3550  0.3880 0.9600878    
## Residuals        40 36.593  0.9148                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 2 :
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## Time_point        4 57.946 14.4864 61.7145 <2e-16 ***
## group             3  1.107  0.3689  1.5715 0.2113    
## Time_point:group 12  4.141  0.3451  1.4703 0.1762    
## Residuals        40  9.389  0.2347                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#interaction plot

#################
#ggplot2 protein
data %>% 
  group_by(Time_point, group) %>% 
  summarise(ave_protein = mean(protein)) -> data2


data2 %>% 
  ggplot() +
  aes(x = Time_point, y = ave_protein, color = group) +
  geom_line(aes(group = group)) +
  geom_point()

#######################
#ggplot2 fat
data %>% 
  group_by(Time_point, group) %>% 
  summarise(ave_fat = mean(fat)) -> data2


data2 %>% 
  ggplot() +
  aes(x = Time_point, y = ave_fat, color = group) +
  geom_line(aes(group = group)) +
  geom_point()

###########
data %>% 
  group_by(group, Time_point) %>% 
  summarise(ave_protein = mean(protein)) -> data2


data2 %>% 
  ggplot() +
  aes(x = group, y = ave_protein, color = Time_point) +
  geom_line(aes(group=Time_point)) +
  geom_point()

###########
data %>% 
  group_by(group, Time_point) %>% 
  summarise(ave_fat = mean(fat)) -> data2


data2 %>% 
  ggplot() +
  aes(x = group, y = ave_fat, color = Time_point) +
  geom_line(aes(group=Time_point)) +
  geom_point()

#boxplot protein 
data %>% 
  group_by(group, Time_point) %>% 
  summarise(ave_protein = mean(protein)) -> data2
data2 %>% ggplot() + aes(x=Time_point,y=ave_protein) + geom_boxplot()

data2 %>% ggplot() + aes(x=group,y=ave_protein) + geom_boxplot()

#boxplot fat
data %>% 
  group_by(group, Time_point) %>% 
  summarise(ave_fat = mean(fat)) -> data2
data2 %>% ggplot() + aes(x=Time_point,y=ave_fat) + geom_boxplot()

data2 %>% ggplot() + aes(x=group,y=ave_fat) + geom_boxplot()

#protein time point

# bar plot

data %>% 
  group_by(group, Time_point) %>% 
  summarise(ave_protein = mean(protein)) -> data2

p <- data2 %>% ggplot() + aes(x=group,y=ave_protein ,fill =group)+
    geom_bar(stat = "identity",position="dodge")+scale_fill_viridis(discrete=TRUE,option = "E") +
    ggtitle("protein Time Point") +
    theme_ipsum() + 
    theme(
      panel.spacing = unit(.1, "lines"),
      strip.text.x = element_text(size = 8)
    ) +
    xlab("") +
    ylab("protein (%)") +
    facet_wrap(~Time_point)

p <- p + theme(axis.text.x = element_text( size = 5,face="bold"))

p

#fat time point

# bar plot

data %>% 
  group_by(group, Time_point) %>% 
  summarise(ave_fat = mean(fat)) -> data2

p <- data2 %>% ggplot() + aes(x=group,y=ave_fat ,fill =group)+
    geom_bar(stat = "identity",position="dodge")+scale_fill_viridis(discrete=TRUE,option = "E") +
    ggtitle("Fat Time Point") +
    theme_ipsum() + 
    theme(
      panel.spacing = unit(.1, "lines"),
      strip.text.x = element_text(size = 8)
    ) +
    xlab("") +
    ylab("fat (%)") +
    facet_wrap(~Time_point)

p <- p + theme(axis.text.x = element_text( size = 5,face="bold"))

p

Analysis of variance

Soroush

2/6/2020