Plotting multiple covariates within MRDS detection functions

[Melissa Costagliola - Ray]

Dear distance sampling group,

I have been using MRDS models (assuming point independence) for point data and have added additional covariates within both the MR and DS parts of the model. I have now compared AICs of all of these models and have a "top model" which has glare and tidal state (ebb / flood) within the MR model and sea state and glare within the DS model.

I am now trying to add these covariate levels to detection function plots as separate lines using add.df.covar.line (i.e. sea state levels, 0-5 and glare levels, 0-3). However, I am struggling with this and was wondering if someone would be able to help me with this. This is my R code:

TopModelFulmar<- ddf(method='io',
                   mrmodel=~glm(link='logit', formula=~distance+GLARE+TidalPhase),
                   dsmodel = ~mcds(key="hr",formula=~distance+SS+GLARE),
                   data=detections, meta.data=list(point=TRUE, convert.units=conversion.factor,adjustment=NULL, left=150))

plot( TopModelFulmar, which=3, showpoints=FALSE) 
add_df_covar_line(TopModelFulmar,data.frame(GLARE=na.omit(unique(detections$GLARE))), pdf=TRUE)

I look forward to hearing from you.

Many thanks,

Melissa

[Eric Rexstad]

Melissa

I've not used the covariate plotting functions in conjunction with a double observer analysis, so I did some experimentation; you'll need to check whether this works for your situation.

The challenge is that the object created by ddf​ for a double observer analysis has lots of components and plotting that object can take on many forms.  Your code indicated that your plot​ call requested the pooled detection which=3​, so far so good.

I think you need to slightly alter your call to add_df_covar_line​ to send only a portion of the ddf​ object to the function.  Here is my example, using the book.tee.data​ supplied with the mrds​ package:

library(mrds)

data(book.tee.data)

region <- book.tee.data$book.tee.region

egdata <- book.tee.data$book.tee.dataframe

samples <- book.tee.data$book.tee.samples

obs <- book.tee.data$book.tee.obs

result.io <- ddf(dsmodel=~mcds(key = "hn", formula=~sex), mrmodel=~glm(~distance),

                 data=egdata, method="io", meta.data=list(width=4))

plot(result.io, showpoints=FALSE, which=3) # pooled detection plot

add_df_covar_line(result.io$ds, data.frame(sex=c(0,1)), lty=1,col=c("red", "green"))

                                 

                 

Note the first argument in the call to add_df_covar_line​; I sent the ds​ element of result.io​, resulting in this graph:

You might also need to build the data.frame that you send to add_df_covar_line​ in a different way, e.g. seq(1, max(detections$GLARE))​ possibly.

Let us know the outcome.

---

From: distance...@googlegroups.com <distance...@googlegroups.com> on behalf of Melissa Costagliola - Ray <melissaco...@gmail.com>
Sent: 29 June 2022 18:32
To: distance...@googlegroups.com <distance...@googlegroups.com>
Subject: {Suspected Spam} [distance-sampling] Plotting multiple covariates within MRDS detection functions

 

--
You received this message because you are subscribed to the Google Groups "distance-sampling" group.
To unsubscribe from this group and stop receiving emails from it, send an email to distance-sampl...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/distance-sampling/CA%2BOf%2Bp4GhUQ8fQGLSLm25ShC-_2J%3DufFyktPUtYHqZ%3DWxBNGkQ%40mail.gmail.com.

[Eric Howe]

Good day Melissa and Eric,

When we include only the ds element of the fitted model, the added lines always reach p = 1 at distance 0 (because we're ignoring the MR model). Eric's example works well because in the book.tee.data, the pooled probability of detection is very close to 1, so the plot looks OK because the originally plotted line reaches p = 1.0. When I try it where p(0) from the MR model is closer to 0.85 the inconsistency between the original plot and the added lines is more apparent and so the plot is less useful.

plot(mrds.model, which = 3, showpoints = F, lwd = 1, breaks = seq(0, 2000, 200))
add.df.covar.line(mrds.model$ds, data = data.frame(SUNANGLE = 1), lty = 2)
add.df.covar.line(mrds.model$ds, data = data.frame(SUNANGLE = 2), lty = 3)
legend(1200, 1, c("Average", "No glare", "Glare"), lty = 1:3, bty = "n")

We can correct the plot to show the effect of the covariate on only the ds model by including only the ds element in the call to ds.plot, but this isn't ideal either when p(0) from the fitted MRDS model is < 1.0.

plot(mrds.model$ds, showpoints = F, lwd = 1, breaks = seq(0, 2000, 200))
add.df.covar.line(mrds.model$ds, data = data.frame(SUNANGLE = 1), lty = 2)
add.df.covar.line(mrds.model$ds, data = data.frame(SUNANGLE = 2), lty = 3)
legend(1200, 1, c("Average", "No glare", "Glare"), lty = 1:3, bty = "n")

I'd like to start with a plot of the MRDS model with p(0) < 1, then subtract 1 - Average combined p(0) from the y-values when calling add.df.covar.line(), but I wasn't able to do that because I wasn't able to save or assign the result of calling plot.ds()

> ds.plot = plot(mrds.model$ds)
> ds.plot
NULL

Any suggestions would be welcome. I may try to calculate y-values from the DS submodel manually (using formulae for half-normal and hazard rate functions and covariate coefficients), then subtract 1 - p(0) from the y-values for plotting.

All the best,

Eric H

[Eric Howe]

Good day,

Of course I was wrong to suggest subtracting 1- average combined p(0) from the y-values predicted from the ds model, because  p(y)=p(0)g(y).

I'm now able to plot the effect of covariates of the scale of the detection function where p(0) is < 1 and estimated from the mr submodel.

E.g. For a ddf object called "fit.br.obsXside.vis.sun.ds.hn0.veg.sun":

# get p(0) from mr submodel
p0.hn = summary(fit.br.obsXside.vis.sun.ds.hn0.veg.sun)$mr.summary$average.p0

# get ds model coefficients
summary(fit.br.obsXside.vis.sun.ds.hn0.veg.sun)$ds$coeff
#$key.scale
#              estimate        se
#(Intercept)  7.5347202 0.3944996
#veg         -0.2411155 0.1215937
#SUNANGLE    -0.4387234 0.2432790

# calculate g(y) for specific values of covariates

hn.det.fn <-function(x, sigma) {
  gy.hn = exp(-x^2/(2*sigma^2))
  return(gy.hn)
}

gy.hn.veg1.sun1 = hn.det.fn(x = seq(0, 2000, 10), sigma = exp(7.5347202 + (1 * -0.2411155) + (1 * -0.4387234)))
gy.hn.veg2.sun1 = hn.det.fn(x = seq(0, 2000, 10), sigma = exp(7.5347202 + (2 * -0.2411155) + (1 * -0.4387234)))
gy.hn.veg3.sun1 = hn.det.fn(x = seq(0, 2000, 10), sigma = exp(7.5347202 + (3 * -0.2411155) + (1 * -0.4387234)))

# calculate p(y) from p(0) and g(y)

py.hn.veg1.sun1 = p0.hn * gy.hn.veg1.sun1
py.hn.veg2.sun1 = p0.hn * gy.hn.veg2.sun1
py.hn.veg3.sun1 = p0.hn * gy.hn.veg3.sun1

plot(fit.br.obsXside.vis.sun.ds.hn0.veg.sun, showpoints = F, which = 3, breaks = seq(0, 2000, 100), main = NULL, lwd = 2)
lines(x = seq(0, 2000, 10), y = py.hn.veg1.sun1, col = "lightgreen", lwd = 2)
lines(x = seq(0, 2000, 10), y = py.hn.veg2.sun1, col = "green", lwd = 2)
lines(x = seq(0, 2000, 10), y = py.hn.veg3.sun1, col = "darkgreen", lwd = 2)
legend(1500, 1, c("Average", "veg=1", "veg=2", "veg=3"), lty=1, lwd = 2,
       col=c("black", "lightgreen", "green", "darkgreen"))

- Eric H