0
$\begingroup$

Our research group studies edge effects of disturbances on adjacent relatively undisturbed ecosystems. We evaluate metrics such as light levels that may be higher at forest edges relative to the interior and the cover of vegetation groups. At the moment, we are exploring continuous data analysis techniques to analyze our data and are considering traditional non-linear least squares (NLS) regression with a specified model (e.g. asymptotic regression) as well as generalized additive modeling (GAM). We aim to determine the magnitude (difference between interior and edge) and the extent (how far it reaches) of the edge effect. We expect that the furthest away distances from the edge should have similar values to each other (reaching an asymptote). However, the curve to reach the asymptote is not necessarily a simple asymptotic regression with a single decline or increase, there might be some local maximum or minimum reached before getting to the final asymptote. In more complex cases, an asymptotic regression would thus not be sufficient and GAM might be more helpful.

This question is closely related to this one from 2021 where the OP modeled a relationship with GAM that suggested an initial increase followed by a decrease (based on distance from location) while the reviewer on the paper insisted that it was overfitting the data when a simple asymptotic regression would have been appropriate. The answers suggest that GAM isn’t just “a polynomial basis of sorts”, however that is how it appears to me and I'd appreciate some guidance on how to coerce GAM to be more than a collection of polynomials.

I have run some simulations where I assumed a simple asymptotic regression curve best described the relationship between the metric and the distance from edge (with normal error terms). The asymptotic regression assumes an artificial extent (distance) of edge effects that is fixed close to ~5m when values reach ~2.3. Then, I analysed the simulated data with an nls regression and with GAM. The output from the nls is as expected. However, the gam curve does not flatten out as expected and instead is quite wiggly beyond the point at which no difference are expected which I cannot fix by alterating parameters (e.g.k, bs). Is this a limitation of GAM or a limitation of my skill-level (e.g. is there a different basis function that should be supplied by me)? When I run the GAM restricting the distance to beyond 5m only, it is a flat curve (non-significant), which is expected. But shouldn't the curve the same on each segment of the data regardless of how it is subset? Increasing the sample size helps match the model to the expected curve and reduces the height of the "bumps" but none of the diagnostic tools that I tried indicate that the model with the smaller sample is problematic (the r-sq is high, the p-values are low, and the residuals seem fine).

gam_values<-gam(nls.rnorm~s(dist_m, k=20), method='REML', data=picked.plt)

GAM with sample size set at 8 (a typical small sample when data is hard to collect such as in the 2021 question example). [blue line is the predefined function used for simulations, black line is the model with sample data, grey points are data input from sample data]

enter image description here enter image description here

GAM with sample size set at 100 (a more generous sample size).

gam_values<-gam(nls.rnorm~s(dist_m, k=3), method='REML', data=picked.plt)

enter image description here enter image description here

GAM with sample size set at 8 - another iteration (random sampling from theoretical) with even worse fit enter image description here

Simulations with a dip (or local minimum) before levelling out enter image description here Full script to run simulations here.

#create set of distances (x-values)

#at 1m intervals
distances_rand<-c(0.5, 1.5, 2.5, seq(3.5, 5.5, by=1), seq(6.5, 10.5, by=1),
                  seq(11.5, 20.5, by=1), seq(21.5, 40.5, by=1), seq(41.5, 80.5, by=1))
dist_incre<-c(1, 2, 3, rep(4, 3), rep(5, 5), rep(6, 10), rep(7, 20), rep(8, 40))
pick_set<-data.table(cbind(distances_rand, dist_incre))

#create set of distances from fixed intervals
picked.plt<-data.table(dist_incre=integer(), distances_rand=integer())

for (i in 1:8){
  picked<-pick_set[, .SD[sample(.N, min(1, .N))], by = dist_incre]
  #picked.pause<-cbind(picked, PlotID=i)
  picked.plt<-rbind(picked.plt, picked)
}

picked.plt<-picked.plt%>%rename("dist_m"="distances_rand")

#predict metric from asymptotic regression NLS pre-defined function
y.preds<-2.31668/(1-0.15144*exp(-(0.44544*picked.plt$dist_m)))
picked.plt$nls.pred<-y.preds

#create simulated normal data with sd=0.1
picked.plt$nls.rnorm<-picked.plt$nls.pred+rnorm(n=nrow(picked.plt), sd=0.1)
#picked.plt$nls.rnorm<-picked.plt$nls.pred+rnorm(n=nrow(picked.plt), sd=0.85)

#run NLS and plot results
nls_values<-nls(nls.rnorm~a/(1+b*exp(-(c*dist_m))), data=picked.plt,
                start=list(a=2.3,b=-0.1, c=0.5))
summary(nls_values)

y.preds2<-summary(nls_values)$coef[1,1]/(1+summary(nls_values)$coef[2,1]*exp(-(summary(nls_values)$coef[3,1]*picked.plt$dist_m)))
picked.plt$nls.pred2<-y.preds2

par(mfrow=c(1,2))
plot(nls.rnorm~dist_m, picked.plt, log="x", main="Output from NLS", col="gray")
lines(nls.pred2~dist_m, picked.plt[order(dist_m),], col="black")#simulated data curve
lines(nls.pred~dist_m, picked.plt[order(dist_m),], col="blue")#initial relationship curve
abline(v=5, col="red")

#run GAM analysis and plot results
gam_values<-gam(nls.rnorm~s(dist_m, k=20), method='REML', data=picked.plt)
#gam_values<-gam(nls.rnorm~s(dist_m, k=3), method='REML', data=picked.plt[dist_m>5.5,])
summary(gam_values)

#plot 
plot(gam_values, log="x", main="Output from GAM", ylim=c(-0.5, 0.5), xlim=c(0.5, 80))
points((nls.rnorm-2.4169)~dist_m, picked.plt, col="gray")
lines(nls.pred-2.4169~dist_m, picked.plt[order(dist_m),], col="blue")#initial relationship curve
abline(v=5, col="red")

par(mfrow=c(2,2))
gam.check(gam_values)

##generate random data for edge effects with local minimum before leveling out####
#input data to obtain a GAM relationship used for simulation
light.dat<-data.table(dist_m = c(0, 0, 0.5, 0.5, 1, 1, 2, 2, 
                                 5, 5, 7, 8, 15,  15,
                                  25, 25, 50, 50, 75, 75),
                       light = c(3.47, 2.71, 3.17, 2.45, 2.8, 2.41, 2.71, 2.26, 
                                 2.13, 1.95, 2.2, 1.8, 2.34, 2.06, 
                                 2.65, 2.09, 2.4, 2.54, 2.56, 2.33))

light.bump<-gam(light~s(dist_m, bs='tp', k=8), light.dat, family="gaussian", method='REML')#smoother straight line beyond edge influence
plot(light.bump, residuals=TRUE, cex=1, pch=1)
o<-plot(light.bump, residuals=TRUE, cex=1, pch=1)

#create set of distances (x-values)
#at 1m intervals
distances_rand<-c(0.5, 1.5, 2.5, seq(3.5, 5.5, by=1), seq(6.5, 10.5, by=1),
                  seq(11.5, 20.5, by=1), seq(21.5, 40.5, by=1), seq(41.5, 80.5, by=1))
dist_incre<-c(1, 2, 3, rep(4, 3), rep(5, 5), rep(6, 10), rep(7, 20), rep(8, 40))
pick_set<-data.table(cbind(distances_rand, dist_incre))

#create set of distances from fixed intervals
picked.plt<-data.table(dist_incre=integer(), distances_rand=integer())

for (i in 1:8){
  picked<-pick_set[, .SD[sample(.N, min(1, .N))], by = dist_incre]
  #picked.pause<-cbind(picked, PlotID=i)
  picked.plt<-rbind(picked.plt, picked)
}

picked.plt<-picked.plt%>%rename("dist_m"="distances_rand")

picked.plt$gam.pred<-predict(light.bump, picked.plt[,"dist_m"])
picked.plt$gam.rnorm<-picked.plt$gam.pred+rnorm(n=nrow(picked.plt), sd=0.2)

#plot(gam.rnorm~dist_m, picked.plt, log="x")

light.bump2<-gam(gam.rnorm~s(dist_m, bs='tp', k=20), picked.plt, 
                 family="gaussian", method='REML')#smoother straight line beyond edge influence
#light.bump2<-scam(gam.rnorm~s(dist_m, bs='mpd', k=20), picked.plt, 
#                 family="gaussian")#smoother straight line beyond edge influence

plot(light.bump2, residuals=TRUE, cex=0.75, pch=1, log="x", shade=TRUE,
     main="GAM, n=8")
     #main="GAM, n=50")
lines(o[[1]]$fit~o[[1]]$x, col="blue")

summary(light.bump2)
$\endgroup$

1 Answer 1

2
$\begingroup$

It might be worth looking into the scam package (https://cran.r-project.org/web/packages/scam/index.html) or something similar, which can set GAMs using monotonically increasing or monotonically decreasing smooth functions. These are still capable of capturing some of the possible complexities that you describe (local humps as you approach the distance threshold), but I would expect them to do a better job of learning when to flatten out.

[example plots added by OP] Run with the simulations from the questions, the output flattens out as expected instead of having multiple humps past the distance of influence.

scam(nls.rnorm~s(dist_m, k=20, bs='mpd'), data=picked.plt)

Output from n=100. Blue line is the theoretical relationship, black line is the modelled relationship. enter image description here

Output from n=8 enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks for your suggestion. On the original hypothetical relationship, it does indeed lead to a smoother outcome -- as I posted the graphs into your answer. Could you provide more suggestions on which options to select if there is a local bump before it flattens out? I am quite bewildered by all the options under shape.constrained.smooth.terms in the package documentation and cannot visualize what the mathematical descriptions should look like. $\endgroup$ Commented Mar 15 at 15:39
  • $\begingroup$ Unfortunately I haven't used these myself so I am not familiar with all the options. But at least those plots suggest that the mpd basis is giving you something useful $\endgroup$ Commented Mar 17 at 23:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.