# downloading data file from Dryad
dryadFileDownload <- function(filenum,filename,baseurl="https://datadryad.org/api/v2")
{
    download.file(paste(baseurl,"/files/",filenum,"/download",sep=""),filename,mode="wb")
}
# load tidyverse
library(tidyverse)
# read in the dataset and show the top few lines
vernFileID <- "22636"
tmpfile <- tempfile(fileext="txt")
dryadFileDownload(vernFileID,tmpfile)
vern <- read.table(tmpfile)
head(vern)

# select the three columns of interest: genotype, vernalization, and phenotype
hua <- select(vern,Genotype,Vernalization,Cauline.leaf.num)
summary(hua)

# filter only observations with `hua2-3` and `hua2-3 FRI` genotypes.
hua <- filter(hua, (Genotype == "hua2-3") | (Genotype == "hua2-3 FRI"))
summary(hua)

# update factor levels
hua <- droplevels(hua)
summary(hua)

# calculate means grouping by genotype and vernalization
huameans <- group_by(hua, Genotype, Vernalization) %>% 
        summarise(meancauline = mean(Cauline.leaf.num, na.rm = T), 
                  n = sum(!is.na(Cauline.leaf.num)))
huameans

6.294118 - 6.453125 # vern effect in FRI

4.761194 - 5.720588 # vern effect in hua2-3

-0.159007 - (-0.959394000000001)

6.453125 - 5.720588	

# group the previous dataframe by genotype
# calculate a weighted mean since the two vernalization groups
# have difference sample sizes
huagenomeans <- group_by(huameans,Genotype) %>% 
                summarize(meancauline = weighted.mean(meancauline, w = n ))
huagenomeans

# the main effect is the difference
huagenomeans$meancauline |> diff() |> round(4)

# the main effect calulated from the hua data
group_by(hua,Genotype) %>% 
    summarize(meancauline = mean(Cauline.leaf.num,na.rm=T))

# interaction term is the difference in the vernalization main effects stratified by genotype
# main effect of vernalization in FRI
( (huameans$meancauline[4] - huameans$meancauline[3]) - 
# main effect of vernalization in hua2-3
  (huameans$meancauline[2] - huameans$meancauline[1]) ) |> 
  round(4) # round to 4 digits

# main effect of genotype
lm(Cauline.leaf.num ~ Genotype, data = hua) |> summary()

# main effect of vernalization
lm(Cauline.leaf.num ~ Vernalization, data = hua) |> summary()

# interaction effect of genotype and vernalization
lm(Cauline.leaf.num ~ Genotype * Vernalization, data = hua) |> summary()

library(lmabc)

lmabc(Cauline.leaf.num ~ Genotype, data = hua) |> summary()

lmabc(Cauline.leaf.num ~ Vernalization, data = hua) |> summary()

lmabc(Cauline.leaf.num ~ Genotype*Vernalization, data = hua) |> summary()

# filter by the six genotypes
vern1 <- filter(vern,(Genotype == "hua2-3") | (Genotype == "hua2-3 FRI")|
                     (Genotype == "vin3-4") | (Genotype == "vin3-4 FRI")|
                     (Genotype == "Col Ama")| (Genotype == "Col FRI")) %>% droplevels()
vern1$Genotype <- factor(vern1$Genotype)
# get counts by genotype
group_by(vern1,Genotype) %>% summarise(n=n())

# filter by not vernalized
# group by genotype
# get mean cauline leaf number (make sure that NA's are removed)
cmeans <- filter(vern1, Vernalization == "NV") %>% 
          group_by(Genotype) %>% 
          summarise(meancauline = mean(Cauline.leaf.num,na.rm=T),
                    n = sum(!is.na(Cauline.leaf.num)))
cmeans

# main effect with the sum contrast is the genotype mean minus mean of all genotype means
dfMainEffect <- mutate(cmeans, maineffect = meancauline-mean(meancauline))
dfMainEffect

levels(as.factor(vern1$Genotype))

# the sum of the main effects is zero
mutate(cmeans,maineffect = meancauline - mean(meancauline)) %>% summarise(sum = sum(maineffect)) %>% round(5) 

options()$contrasts

# change contrast option
options(contrasts=c("contr.sum","contr.poly"))
# use lm to get main effects
(out1 <- lm(Cauline.leaf.num ~ Genotype, data = vern1[vern1$Vernalization=="NV",])) %>% summary

lmabc(Cauline.leaf.num ~ Genotype, data = vern1[vern1$Vernalization=="NV",]) |> summary()

options()$contrasts

# all but one main effect
(meff <- t(coef(out1)[-1])) %>% round(4)

# negative of the sum
-sum(meff) %>% round(4)

# relevel
vern1$Genotype <- relevel(vern1$Genotype,ref="vin3-4 FRI")

# fit new model
out2 <- lm(Cauline.leaf.num ~ Genotype, 
           data = vern1[vern1$Vernalization=="NV",])
summary(out2)

data.frame(original=coef(out1),releveled=coef(out2))