FirmnessProject_Coding_Final_Adjusted.R

install.packages("dplyr")
install.packages("cowplot")
install.packages("psych")
library(cowplot)
library(dplyr)
library(car)
library(agricolae)
library(ExpDes)
library(ggplot2)
library(psych)

firm <- read.csv("C:/Users/RhysB/OneDrive/Desktop/Texture/Rhys Discovery/RSPREAD.csv")
mutate(firm, harv_no = factor(harv_no), rep = factor(rep))
firm

#ANOVA will be performed to determine differences in fruit firmness among genotypes
#for each method with harvest dates as replicates.

#Filter data by type

ff_subset <- subset(firm, type == "ff")
taxt_subset <- subset(firm, type == "taxt")
subj_subset <- subset(firm, type == "subj")

#ANOVA for each type/method

ff.firm.aov <- aov(force ~ seln, data = ff_subset)
summary(ff.firm.aov)

taxt.firm.aov <- aov(force ~ seln, data = taxt_subset)
summary(taxt.firm.aov)

subj.firm.aov <- aov(force ~ seln, data = subj_subset)
summary(subj.firm.aov)

time.aov <- aov(time ~ type, data = firm)
summary(time.aov)
#TIME CHART
# Calculate average time for each method
avg_time <- firm %>%
  group_by(type) %>%
  summarise(avg_time = mean(time))

# Create a bar chart

ggplot(avg_time, aes(x = reorder(type, avg_time), y = avg_time)) +
  geom_bar(stat = "identity", fill = "black", color = "black", size = 0.05) +
  labs(x = "Test Type", y = "Mean Time (seconds/berry)") +
  scale_x_discrete(labels = c("ff", "subj", "taxt")) +
  coord_flip() +
  theme_minimal() +
  theme(panel.grid.major = element_line(size = 0.5, color = "#F7F7F7"),
        panel.grid.minor = element_line(size = 0.5, color = "#F7F7F7"),
        axis.line = element_line(size = 0.5, color = "#455A64"),
        axis.text = element_text(size = 12, color = "#455A64"),
        axis.title = element_text(size = 14, color = "black"),
        axis.title.y = element_text(margin = margin(r = 10)))

#All have significant p-values, aka all models find differences in firmness values between genotypes
#Post-Hoc Tukey Analysis for each ANOVA model
#The error values come from the associated ANOVA output above
#Mean separation will be performed with Tukey’s Honestly Significant Difference where appropriate.




time.hsd <- HSD.test(time.aov, "type", group=TRUE)
print(time.hsd)

tukey(ff_subset$force, ff_subset$seln, DFerror=470, SSerror=358557)

tukey(taxt_subset$force, taxt_subset$seln, DFerror=470, SSerror=25272682)

tukey(subj_subset$force, subj_subset$seln, DFerror=470, SSerror=378.8)

#check out the order between these 3, this is what we will make the box plots with
#the order is similar between ff and taxt, but subj is all over the place
#One-way ANOVA will be performed to detect significant differences in time required for
#taxt takes way longer on average than subj or ff
#This time difference along with the better correlation between ff and taxt is your justification for using ff in the future
#ff is better related to taxt than subj is
#ff is faster than taxt

#ordering from low to high firmness values, you may want to keep this order for each figure so you can better visualize the variability
#ff_subset$seln <- reorder(ff_subset$seln, ff_subset$force, FUN = median)
#don't do this because it will change the order in the boxplots
#this code creates the plot, so you can copy this and modify it for each additional figure

a <- ggplot(ff_subset, aes(x = seln, y = force)) +
    geom_boxplot()  + labs(x = "Genotype", y = "Force") + ggtitle("FruitFirm 1000") +
  geom_text(data = ff_subset, aes(label = meansep, y = ff_subset$force, vjust=-0.3), 
            position = position_dodge(width = 0.75), 
            show.legend = FALSE )
a

#TAXT boxplot
b <- ggplot(taxt_subset, aes(x = seln, y = force)) +  geom_boxplot() +
    labs(x = "Genotype", y = "Force") +  ggtitle("TA.XT") +
  geom_text(data = taxt_subset, aes(label = meansep, y = taxt_subset$force, vjust=-0.3), 
            position = position_dodge(width = 0.75), 
            show.legend = FALSE )
b

#Subj boxplot
c <- ggplot(subj_subset, aes(x = seln, y = force)) +  geom_boxplot() +
    labs(x = "Genotype", y = "Force") +  ggtitle("Subjective") +
  geom_text(data = subj_subset, aes(label = meansep, y = subj_subset$force, vjust=-0.3), 
            position = position_dodge(width = 0.75), 
            show.legend = FALSE ) + expand_limits(y=6)
c


#Need something here to widen row barriers so we can see the whole mean sep letter
a <- ggplot(ff_subset, aes(x = seln, y = force)) +
  geom_boxplot()  + labs(x = "Genotype", y = "Force") + ggtitle("FruitFirm 1000") +
  geom_text(data = ff_subset, aes(label = meansep, y = ff_subset$force + 0.12*max(ff_subset$force)), 
            position = position_dodge(width = 0.75), 
            show.legend = FALSE ) +
  theme(panel.spacing.y = unit(1, "lines"))

b <- ggplot(taxt_subset, aes(x = seln, y = force)) +  geom_boxplot() +
  labs(x = "Genotype", y = "Force") +  ggtitle("TA.XT") +
  geom_text(data = taxt_subset, aes(label = meansep, y = taxt_subset$force + 0.12*max(taxt_subset$force)), 
            position = position_dodge(width = 0.75), 
            show.legend = FALSE ) +
  theme(panel.spacing.y = unit(1, "lines"))

c <- ggplot(subj_subset, aes(x = seln, y = force)) +  geom_boxplot() +
  labs(x = "Genotype", y = "Force") +  ggtitle("Subjective") +
  geom_text(data = subj_subset, aes(label = meansep, y = subj_subset$force + 0.12*max(subj_subset$force)), 
            position = position_dodge(width = 0.75), 
            show.legend = FALSE ) +
  theme(panel.spacing.y = unit(1, "lines")) +
  expand_limits(y=6)

#using cowplot to merge 3 plots into one

plot_grid(a, b, c, labels = c("a.", "b.", "c."), ncol = 1, nrow = 3, 
          rel_heights = c(1, 1, 1), 
          heights = unit(c(1, 1, 1), "null", 1.5))

#Pearson’s correlations for firmness measurements between methods will be tested using PROC CORR (SAS 9.4, Cary, NC).



# Calculate average firmness per genotype and measurement type

average_firmness <- aggregate(force ~ seln + type, data = firm, FUN = mean)



# Reshape the data to have each measurement type as a separate column

average_firmness_wide <- reshape(average_firmness, idvar = "seln", timevar = "type", direction = "wide")
average_firmness_wide

firm

firm_wide <- reshape(firm, idvar = "seln", timevar="type", direction="wide")

# Compute correlations between average firmness for different measurement types

correlations <- cor(average_firmness_wide[, 2:ncol(average_firmness_wide)], method = "pearson")
correlations

# Print correlation matrix

print(correlations)

pairs(average_firmness_wide[,2:4], pch = 19)

pairs(average_firmness_wide[,2:4], pch = 19,  cex = 1.0,
      lower.panel=NULL)

panel.cor <- function(x, y){
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- round(cor(x, y), digits=2)
  txt <- paste0("R = ", r)
  text(0.5, 0.5, txt, cex = 1.5*par("cex.lab"), font = 1)  # Increase font size of R= boxes
}  


upper.panel<-function(x, y){
  points(x,y, pch = 19)
}

pairs(average_firmness_wide[,2:4], 
      lower.panel = panel.cor,
      upper.panel = upper.panel)

)