In mdlama/milkweed:

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

library(milkweed)
library(dplyr)
library(ggplot2)
library(fitdistrplus)
`%$%` <- magrittr::`%$%`

Read, filter and transform data.

metadata <- tbl_df(stemdata)

Fits of h_apical for seedlings.

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

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)

Looks like Weibull is the best fit, but similar to budling distributions, we'll go with normal for easier interpretation with perturbation analysis.

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