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)
library(magrittr)
sites <- c("Bertha", levels(stemdata$site))

Read and filter.

metadata <- tbl_df(stemdata)

metadata %<>% filter(!is.na(h_apical),
                     !is.na(munched))

Bertha

First grab probability of being munched.

(pmunch <- sum(metadata$munched == 1)/nrow(metadata))

Now do fits for herbivory scores of munched plants.

herb_avg <- (metadata %>% filter(munched == 1))$herb_avg

f0 <- fitdist(herb_avg, "norm")
f1 <- fitdist(herb_avg, "lnorm")
f2 <- fitdist(herb_avg, "gamma")
f3 <- fitdist(herb_avg, "exp")
f4 <- fitdist(herb_avg, "weibull")

Plot results.

par(mfrow=c(2,2))
plot.legend <- c("norm", "lnorm", "gamma", "exp", "weibull")
denscomp(list(f0, f1, f2, f3, f4), legendtext = plot.legend)
cdfcomp(list(f0, f1, f2, f3, f4), legendtext = plot.legend)
qqcomp(list(f0, f1, f2, f3, f4), legendtext = plot.legend)
ppcomp(list(f0, f1, f2, f3, f4), legendtext = plot.legend)

Goodness-of-fit measures.

gofstat(list(f0, f1, f2, f3, f4), fitnames = plot.legend)

Estimates. The inputs to predict always on [0, xmax], broken up into equal intervals.

munched.fit <- vector('list', length(sites))
names(munched.fit) <- sites

munched.fit[['Bertha']] <- vector("list", 3)
munched.fit[['Bertha']][[1]] <- f1
munched.fit[['Bertha']][[2]] <- pmunch

munched.fit[['Bertha']][[3]] <- eval(parse(text = sprintf("function(x, justmunch=FALSE) {
  N <- length(x)
  dx <- x[2]-x[1]
  y <- rep(0, N-1)
  for (j in 1:(N-1)) {
    y[j] = p%s(x[j+1], %g, %g) - p%s(x[j], %g, %g)
  }
  y[1] <- y[1] + p%s(x[1], %g, %g)
  y[N-1] <- y[N-1] + p%s(x[N], %g, %g, lower.tail=FALSE)
  y <- %g*y
  if (!justmunch) {
    y[1] <- y[1] + (1-%g)
  }
  y <- y/dx
}",
  munched.fit[[1]][[1]]$distname, munched.fit[[1]][[1]]$estimate[1], munched.fit[[1]][[1]]$estimate[2],
  munched.fit[[1]][[1]]$distname, munched.fit[[1]][[1]]$estimate[1], munched.fit[[1]][[1]]$estimate[2],
  munched.fit[[1]][[1]]$distname, munched.fit[[1]][[1]]$estimate[1], munched.fit[[1]][[1]]$estimate[2],
  munched.fit[[1]][[1]]$distname, munched.fit[[1]][[1]]$estimate[1], munched.fit[[1]][[1]]$estimate[2],
  munched.fit[[1]][[2]], munched.fit[[1]][[2]]
)))
names(munched.fit[['Bertha']]) <- c("fit", "pmunch", "predict")

To illustrate how this would be used, we need a point mass at the origin with the rest the log-normal.

par(mfrow=c(1,1))
bx <- seq(from = 0, to = 6, length.out=101)
x <- 0.5*(bx[1:100] + bx[2:101])
herb_avg <- metadata$herb_avg
hist(herb_avg, breaks=seq(from = 0, to = 6, by=0.05), freq=FALSE)
lines(x, munched.fit[['Bertha']]$predict(bx), col="red", lwd=2)

Sites

# munched.fit[["PWR"]]$fit$estimate <- c(-2.25, 1.1)

for (i in 2:length(sites)) {
  thissite <- metadata %>% filter(site == sites[i])
  pmunch <- sum(thissite$munched == 1)/nrow(thissite)
  herb_avg <- (thissite %>% filter(munched == 1))$herb_avg

  f0 <- fitdist(herb_avg, "norm")
  f1 <- fitdist(herb_avg, "lnorm")
  f2 <- fitdist(herb_avg, "gamma")
  f3 <- fitdist(herb_avg, "exp")
  f4 <- fitdist(herb_avg, "weibull")

  # ind <- which.min(gofstat(list(f0, f1, f2, f3, f4))$aic)
  ind <- 2 # We'll use log-normal for all

  munched.fit[[i]] <- vector("list", 3)
  eval(parse(text=sprintf("munched.fit[[%d]][[1]] <- f%d", i, ind-1)))
  # munched.fit[["PWR"]][[1]]$estimate <- c(-2.25, 1.1)
  munched.fit[[i]][[2]] <- pmunch
  munched.fit[[i]][[3]] <- eval(parse(text = sprintf("function(x, justmunch=FALSE) {
  N <- length(x)
  dx <- x[2]-x[1]
  y <- rep(0, N-1)
  for (j in 1:(N-1)) {
    y[j] = p%s(x[j+1], %g, %g) - p%s(x[j], %g, %g)
  }
  y[1] <- y[1] + p%s(x[1], %g, %g)
  y[N-1] <- y[N-1] + p%s(x[N], %g, %g, lower.tail=FALSE)
  y <- %g*y
  if (!justmunch) {
    y[1] <- y[1] + (1-%g)
  }
  y <- y/dx
}",
  munched.fit[[i]][[1]]$distname, munched.fit[[i]][[1]]$estimate[1], munched.fit[[i]][[1]]$estimate[2],
  munched.fit[[i]][[1]]$distname, munched.fit[[i]][[1]]$estimate[1], munched.fit[[i]][[1]]$estimate[2],
  munched.fit[[i]][[1]]$distname, munched.fit[[i]][[1]]$estimate[1], munched.fit[[i]][[1]]$estimate[2],
  munched.fit[[i]][[1]]$distname, munched.fit[[i]][[1]]$estimate[1], munched.fit[[i]][[1]]$estimate[2],
  munched.fit[[i]][[2]], munched.fit[[i]][[2]]
  )))

  names(munched.fit[[i]]) <- c("fit", "pmunch", "predict")
  cat("\n-------\n")
  cat(sites[i], "\n")
  print(munched.fit[[sites[i]]])
}

N = 50
lmin = 0
lmax = 6
b = lmin + c(0:N)*(lmax - lmin)/N # boundary points
x = 0.5*(b[1:N] + b[2:(N+1)])
dx = b[2] - b[1] # class size

curves <- do.call(rbind, lapply(sites[-1],
                                function(site) {
                                  data.frame(site = site,
                                             x = x,
                                             y = munched.fit[[site]]$predict(b))
                                }
))

p <- metadata %>% ggplot(aes(x = herb_avg)) +
  geom_histogram(aes(y = ..density..),
                 binwidth = 0.1) +
  geom_line(aes(x = x, 
                y = y), 
            data = curves, 
            colour = "red") +
  facet_wrap(~ site) +
  coord_cartesian(ylim=c(0, 3)) +
  theme_bw() +
  xlab("Herbivory Score") +
  ylab("Density")

p

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