I have to fit three regressions (an ordinary linear model, a beta-regression (offered by betareg’ betareg) and another beta-regression (offered by mgcv’s family option betarto)) with the data of the cover of blueberry in pine forest in northern Sweden which is somewhat affected by the thickness
of the soil humus layer. Both variables are in two different data sets.
Also I need to plot the three fitted models, including 95% confidence interval, into a “production-level” plot.
This is what I have. Can somebody help me with the code?
Code: Alles auswählen
install.packages("vegan")
install.packages("betareg")
install.packages("mgcv")
library(vegan)
data(varechem)
data(varespec)
varespec$NewBlueberryCover <- varespec$Vaccmyrt + 0.1
varechem$Humdepth <- varechem$Humdepth / 100
varespec$NewBlueberryCover <- varespec$NewBlueberryCover / 100
plot(varechem$Humdepth, varespec$NewBlueberryCover,
xlab = "Humus Depth", ylab = "Blueberry Cover",
main = "Blueberry Cover with Humus Depth")
combined_data <- data.frame(Humdepth = varechem$Humdepth, NewBlueberryCover = varespec$NewBlueberryCover)
lm_model <- lm(NewBlueberryCover ~ Humdepth, data = combined_data)
summary(lm_model)
beta_model <- betareg(NewBlueberryCover ~ Humdepth, data = combined_data)
summary(beta_model)
beta_model_mgcv <- gam(NewBlueberryCover ~ s(Humdepth), data = combined_data, family = betar(link = "logit"))
summary(beta_model_mgcv)
cols <- terrain.colors(4)
col2rgb(cols)
plot(NewBlueberryCover ~ Humdepth, data=combined_data, pch=16, col=cols, las=1, xlab="Humus Depth", ylab="Blueberry Cover")
with(combined_data, tapply(NewBlueberryCover, Humdepth))
preds1 <- predict(lm_model, newdata=data.frame("Humdepth"=seq(min(combined_data$Humdepth), se.fit=T) +
lines(seq(min(combined_data$Humdepth), exp(preds1$fitted.values[,1]), col=cols[1], lwd=2, lty=2) +
polygon(c(Humdepth, rev(Humdepth)), c(exp(preds1$fitted.values[,1] + 2*preds1$se.fit[,1]),
rev(exp(preds1$fitted.values[,1] -2*+ 2*preds1$se.fit[,1]))), border=NA, col=rgb(0, 166, 0, 30, maxColorValue=255)))))--> Also I am not sure to combine the variables from the two data sets into one.