R/Yield_functions.R

In henganl2/Yield: Soybean Yield Data Analysis

#Get yield (bushels/acre) data

adjusted_weight <- function(weight, moisture, standard_moisture, standard_bushel_weight, ...){
  # the weight(lb) and moisture are collexted by combine 
  # this is aim to adjust the weight by moisture
  
  if (length(which(is.na(weight)==TRUE)> 0)){
    print("Missing value found in the weight column")
  }
  if (length(which(is.na(moisture)==TRUE)> 0)){
    print("Missing value found in the moisture column")
  } 
  if(any(weight <= 0 | moisture<=0, na.rm=TRUE)){
    print("Unreasonable number of weight or moisture found")
  }
  # A: adjust weight by moisture
  A <-((100-moisture)*0.01)*weight
  # B: standard bushel wieight of soybean 
  B <- ((100-standard_moisture)*0.01)*standard_bushel_weight
  #C: plot_width(ft)*plot_length(ft)/43560
  C <- 7*15/43560
  return(A/B/C)
}

#test normality assumption
ANOVA_norm_assumption <- function(adjusted_weight, Treatment,...){
  ANOVA.model <- lm(adjusted_weight ~ as.factor(Treatment))
  output <- anova(ANOVA.model)
  residuals <- ANOVA.model$residuals 
  normality <-shapiro.test(residuals)
  pvalue<-normality$p
  if(pvalue>0.05){
    cat("data are normally distributed")
  }else{
    cat("data are not normally distributed")
  }
  QQplot <-qqnorm(residuals)
  print(pvalue)
  write.csv(output, file="ANOVA_table.csv")
}


#test HOV assumption 
ANOVA_HOV_assumption <- function(adjusted_weight, Treatment,...){
  ANOVA.model <- lm(adjusted_weight ~ as.factor(Treatment))
  output <- anova(ANOVA.model)
  residuals <- ANOVA.model$residuals
  HOVtest <- levene.test(adjusted_weight, Treatment)
  pvalue<-HOVtest$p
  if(pvalue>0.05){
    cat("HOW are met")
  }else{
    cat("HOV not met ")  
  }
  print(pvalue)
  predicted.values <- ANOVA.model$fitted.values
  plot(residuals~predicted.values)
  abline(h = 0, lty = 2)
  write.csv(output, file="ANOVA_table.csv")
}



#generate some figures
Yield_plot <- function(adjusted_weight,Treatment,...){
par(mfrow=c(1,2))
hist(adjusted_weight, main="Histogram of adjusted_weight",
     xlab="Yield(bushels/acre")
boxplot(adjusted_weight~Treatment,
        xlab="Treatment",
        ylab="yield (bushels/acre)",
        main="Boxplot of adjusted_weight")
}