In mdlama/milkweed:

We use the fdistrplus package (https://cran.r-project.org/web/packages/fitdistrplus/vignettes/paper2JSS.pdf).

Summary While many of the fits suggest going with the Weibull distribution, we will use the normal distribution throughout as the fits are good with the normal distribution, and for ease of implementation and interpretation.

library(milkweed)
library(dplyr)
library(ggplot2)
library(fitdistrplus)
`%$%` <- magrittr::`%$%`
sites <- c("Bertha", levels(stemdata$site))

Read, filter and transform data.

metadata <- tbl_df(stemdata)

h_apical <- metadata %>% 
  filter(seedling == 0,
         !is.na(h_apical)) %$%
  h_apical

Bertha

Fits of h_apical for budlings.

distnames <- c("norm", "lnorm", "gamma", "exp", "weibull")

f0 <- fitdist(h_apical, distnames[1])
f1 <- fitdist(h_apical, distnames[2])
f2 <- fitdist(h_apical, distnames[3])
f3 <- fitdist(h_apical, distnames[4])
f4 <- fitdist(h_apical, distnames[5])

Plot results.

plot_diagnostics <- function(...) {
  par(mfrow=c(2,2))
  denscomp(list(...), legendtext = distnames)
  cdfcomp(list(...), legendtext = distnames)
  qqcomp(list(...), legendtext = distnames)
  ppcomp(list(...), legendtext = distnames)
  par(mfrow=c(1,1))
}

plot_diagnostics(f0, f1, f2, f3, f4)

Goodness-of-fit measures.

gofstat(list(f0, f1, f2, f3, f4), fitnames = distnames)

Weibull appears to beat out normal distribution, but based on the diagnostic plots, normal is good enough as far as we're concerned - makes for easier interpretation and implementation of parameter perturbation for elasticity analysis.

Sites

fit_site <- function(mysite) {
  cat(mysite)
  h_apical <- metadata %>% 
    filter(seedling == 0,
           !is.na(h_apical),
           site == mysite) %$%
    h_apical

  f0 <- fitdist(h_apical, distnames[1])
  f1 <- fitdist(h_apical, distnames[2])
  f2 <- fitdist(h_apical, distnames[3])
  f3 <- fitdist(h_apical, distnames[4])
  f4 <- fitdist(h_apical, distnames[5])

  cat("\n-------\n")
  print(gofstat(list(f0, f1, f2, f3, f4), fitnames = distnames))
  ind <- which.min(gofstat(list(f0, f1, f2, f3, f4))$bic)

  cat("Diagnostic plots:\n")
  plot_diagnostics(f0, f1, f2, f3, f4)
}

BLD1

fit_site("BLD1")

Results suggest going with a normal distribution.

BLD2

fit_site("BLD2")

Results suggest going with a Weibull or normal distribution.

GET

fit_site("GET")

Results suggest going with a normal distribution.

GRN

fit_site("GRN")

Results suggest going with a normal or Weibull distribution.

PWR

fit_site("PWR")

Results suggest going with a normal distribution.

SKY

fit_site("SKY")

Results suggest going with a Weibull distribution.

YTB

fit_site("YTB")

Results suggest going with a Weibull distribution.

mdlama/milkweed documentation built on Sept. 15, 2022, 9:24 a.m.