The data (dry foliage biomass in kg) are not count data, so you shouldn't use either Poisson or quasi-Poisson responses.
The first thing to try is usually a regular linear model:
m1 <- lm(Foliage ~ DBH+Age+Origin, data = lime)
library(performance)
check_model(m1)
This is pretty terrible (if the first row [nonlinearity and heteroscedasticity] is bad, you don't even really need to look at the rest):
The next thing to try for continuous data are positive (e.g. kg of foliage) is to log-transform the response variable:
m2 <- update(m1, log(Foliage) ~ .)
check_model(m2)
This looks much better. There still seems to be slight evidence of nonlinearity and heteroscedasticity, but it's not too bad.
One possible solution is to see if interactions explain any of the patterns. Explore the pattern in the residuals:
library(broom)
aa <- augment(m2)
library(ggplot2)
ggplot(aa, aes(DBH,.std.resid, colour = Origin)) +
geom_point() +
geom_smooth() +
facet_wrap(~cut_number(lime$Age, 4))
The fact that the downturn in the residuals occurs consistently across origins for young/small trees makes me think that it's actually a real pattern. What you would do about this would depend on how badly you needed an accurate model: for example, you could exclude the young/small trees from your model, or fit a generalized additive model, or a nonlinear model ...
As mentioned in the comments and in @GordonSmyth's answer, you could also use a Gamma GLM (probably with a log link), but the results would probably be very similar to the log-linear model ... this is
m3 <- glm(formula = Foliage ~ DBH + Age + Origin,
family = Gamma(link = "log"),
data = lime)
Computing the AIC of the log-linear model (accounting for the transformation of the response variable) and the Gamma model:
AIC(m2) + sum(2*log(lime$Foliage)) - AIC(m3)
## [1] 15.43352
suggests that the Gamma model is indeed better (lower AIC is better; delta-AIC > 10 is a very large difference).
rr <- purrr::map_dfr(list(loglin = m2, Gamma = m3), augment, .id = "model")
ggplot(rr, aes(.fitted, .std.resid, colour = model)) +
geom_point() + geom_smooth()
However, the same nonlinear patterns are still present in the residuals (maybe a little less pronounced for the Gamma fit).
Finally, adding a squared DBH term to the model (m4 <- update(m3, . ~ . + I(DBH^2)
) appears to help a lot ...