p and esw for each value of covariat

[lea david]

Dear,

I use Distance 7.3 software, and I ran the MCDS engine with a covariate. This factor influence the detectability. I would like to know for each of the value of the covariate, the associated p (probabilité of detection) and esw. Where do I find those results in this software ?

Many thanks

Léa

[Eric Rexstad]

Léa

I'm afraid the estimates you suggest are not produced by the Distance for Windows software.  To produce such estimates would require calculation on your part.  You would need to compute the detection function for each level of the covariate from the parameter estimates shown in the Distance for Windows output such as shown in the attached screen grab for a hazard rate detection function that includes a three-level factor covariate.

After computing the "s" or scale parameter for a covariate factor level, then you would substitute that into the "hazard rate key" function.

Next that function would need to be integrated from 0 to your truncation distance and appropriately scaled, to provide an estimate of P_a for that covariate level.

Multiplying the truncation distance "w" by the derived estimate of P_a would deliver an estimate of effective strip width.

If this is a challenge, send me more details off-list and I'll try to offer some help.

---

From: distance...@googlegroups.com <distance...@googlegroups.com> on behalf of lea david <leatti...@gmail.com>
Sent: 11 February 2022 08:34
To: distance-sampling <distance...@googlegroups.com>
Subject: {Suspected Spam} [distance-sampling] p and esw for each value of covariat

 

--
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/e6542c01-508a-45a0-b311-efd70c872d7cn%40googlegroups.com.

[Tim Awbery]

Hi,

I am trying to compute the probability of detection for different sea states in a problem similar to Léa's above (but in R rather than Distance for Windows) in order to obtain the ESW for each sea state.

I understand that the ESW is estimated by the average probability * truncation distance and am able to calculate this for the overall model using the script that Eric helpfully provided on another thread:

ESW.calc <- function(my.ds.object) {
  #  calculates effective strip width
  #    given object created by function ds() of class "dsmodel"
  #    Rexstad, 28Apr17
  if (class(my.ds.object) != "dsmodel") stop("Argument is not of class dsmodel")
  my.summary <- summary(my.ds.object)
  truncation <- my.ds.object$ddf$meta.data$width
  average.p <- my.summary$ds$average.p
  esw <- truncation * average.p
  return(esw)
}

I also understand that when including covariates, you are effectively scaling the shape of the detection function. 

I believe I have computed the scale parameter for each covariate factor level as Eric suggested above:

Distance range : 0 - 733.5 Detection function parameters Scale coefficient(s): estimate se (Intercept) 5.5619879 0.2546760 as.factor(beaufort)1 -0.6515770 0.3216016 as.factor(beaufort)2 -0.6531300 0.2758832 as.factor(beaufort)3 -0.8252728 0.3208008 as.factor(beaufort)4 -1.4570742 0.4011394 Shape coefficient(s): estimate se (Intercept) 0.4260557 0.09593507 Estimate SE CV Average p 0.3242452 0.02433436 0.07504923

What I am struggling with is, is how substitute these values into the hazard rate key function as suggested above and indeed how one would go about integrating the area between 0 and the truncation distance. I have naively tried using the integrate function the "dsmodel" object but obviously this is the whole model rather than the actual function.

I'd be grateful for any advice.

Best,

Tim

[Eric Rexstad]

Tim

A bit of a mess with code to do some fiddly calculations in answer to your questions; borrowing some code from my colleague Tiago along the way.  Need to extract estimated coefficients, place them into the hazard rate key function, integrate said function over the appropriate range to get ESW (according to this formula):

Code has not been checked to ensure its integrity, but it seems to give reasonable results.  Note the coefficients for sea states 1 and 2 are virtually identical, hence the two detection function curves live atop one another in this plot:

gz<-function(z,

             beta, sigintercept, sigcoef, DistWin,

             key="HR", w=max(z)){

  #this is a generic detection function that returns the probability of detecting an animal

  #    z               generic distance (perpendicular) - can be scalar or vector

  #    beta            shape coefficient

  #    sigmaintercept  intercept coefficient for sigma

  #    sigcoef         coefficient for specific factor level

  #    DistWin         coefficients from Distance for Windows or from R

  #    key             the detection function key, only works for hazard rate

  #    w               truncation distance, by default the max of the distances

  #

  #RETURNS: a probability

  if(key=="HR") {

    if (DistWin) {

      sigma <- sigintercept + exp(sigcoef)

      exponent <- beta

    } else {

      sigma <- exp(sigintercept + sigcoef)

      exponent <- exp(beta)

    }

    scale.dist <- z/sigma

    inside <- -(scale.dist)^exponent

    gx <- exp(inside)

  }

  return(gx)

}

beaustates <- c("State 0","State 1", "State 2", "State 3", "State4")

beaucoeff <- c(0, -.6515770, -0.6531300, -0.8252728, -1.4570742)

esw <- vector(mode="numeric", length=length(beaustates))

for (i in 1:length(beaustates)) {

  esw[i] <- integrate(gz, 0, 500,

                      beta=.4260557, sigintercept = 5.5619879, sigcoef = beaucoeff[i],  DistWin=FALSE)$value

  if(i==1) {

    plot(seq(0,500),

         gz(beta=.4260557, sigintercept = 5.5619879, sigcoef = beaucoeff[i], z=seq(0,500), DistWin=FALSE),

         lwd=3, type="l", col=i, main="Det fns by sea state", xlab="Distance", ylab="Detection probability")

  } else {

    lines(seq(0,500),

          gz(beta=.4260557, sigintercept = 5.5619879, sigcoef = beaucoeff[i], z=seq(0,500), DistWin=FALSE),

          lwd=3, type="l", col=i)

  }

}

msg <- paste(beaustates, title="Estimated ESW", round(esw, 1))

legend("topright", legend=msg, lwd=2, col=1:5)

---

From: distance...@googlegroups.com <distance...@googlegroups.com> on behalf of Tim Awbery <tim.a...@sams.ac.uk>
Sent: 10 May 2023 16:42
To: distance-sampling <distance...@googlegroups.com>
Subject: Re: [distance-sampling] p and esw for each value of covariat

 

To view this discussion on the web visit https://groups.google.com/d/msgid/distance-sampling/6b36f004-ab3c-4a72-ba1a-9db1daad03b1n%40googlegroups.com.

[Tim Awbery]

Hi,

 

Thanks for your prompt response and for the provided code.

The function appears to run, and I am able to reproduce the figure above with no problem.

 

By my calculations, my ESW for the overall model is: mu = W * Pa  = 800 * 0.41 = 331.43 metres. 

I would therefore expect my ESW in lower sea states to be greater than this and vice versa in higher sea states.

 However, when I substitute in my model intercept, scale coefficients and shape coefficient, even at a Sea State of zero, my ESW is only 217 metres with the ESW reduced to 58 metres at a Sea State of 5.

 Additionally, the functions appear to have a much smaller shoulder compared to the original model. I thought I understood which values needed replacing (please see attached) but now I am less sure.

 

If you’re able to provide any confirmation or correction regarding this, then I’d be very grateful.

I have attached the output figure provided by yours and Tiago’s code as well as the shape of the original model and my edited code, in case it is of any help in diagnosing my mistake.

 

Thank you in advance.

 

Best,

 

Tim

beaustates <- c("State 0","State 1", "State 2", "State 3", "State4")

beaucoeff <- c(0, -0.3979885, -0.5518057, -0.8339985, -1.3146111)

esw <- vector(mode="numeric", length=length(beaustates))

for (i in 1:length(beaustates)) {

        esw[i] <- integrate(gz, 0, 800,

                            beta=.4260557, sigintercept = 5.4885488, sigcoef = beaucoeff[i],  DistWin=FALSE)$value

        if(i==1) {

                plot(seq(0,800),

                     gz(beta=0.4142978, sigintercept = 5.4885488, sigcoef = beaucoeff[i], z=seq(0,800), DistWin=FALSE),

                     lwd=3, type="l", col=i, main="Det fns by sea state", xlab="Distance", ylab="Detection probability")

        } else {

                lines(seq(0,800),

                      gz(beta=.4260557, sigintercept = 5.4885488, sigcoef = beaucoeff[i], z=seq(0,800), DistWin=FALSE),

                      lwd=3, type="l", col=i)

        }

}

msg <- paste(beaustates, title="Estimated ESW", round(esw, 1))

legend("topright", legend=msg, lwd=2, col=1:5)

[Eric Rexstad]

Tim

Thanks for double-checking the code. Indeed the code did have a gremlin, a messy bunch of minus signs wandering about inside the hazard rate key function. I think this code produces a result that passes the sense check.

gz<-function(z,

             beta, sigintercept, sigcoef, DistWin,

             key="HR", w=max(z)){

  #this is a generic detection function that returns the probability of detecting an animal

  #    z               generic distance (perpendicular) - can be scalar or vector

  #    beta            shape coefficient

  #    sigmaintercept  intercept coefficient for sigma

  #    sigcoef         coefficient for specific factor level

  #    DistWin         coefficients from Distance for Windows or from R

  #    key             the detection function key, only works for hazard rate

  #    w               truncation distance, by default the max of the distances

  #

  #RETURNS: a probability

  if(key=="HR") {

    if (DistWin) {

      sigma <- sigintercept + exp(sigcoef)

      exponent <- beta

    } else {

      sigma <- exp(sigintercept + sigcoef)

      exponent <- exp(beta)

    }

    scale.dist <- z/sigma

    inside <- -(scale.dist)^(-exponent)

    gx <- 1 - exp(inside)

  }

  return(gx)

}

maxz <- 800

betaval <- 0.4142978

sigint <- 5.4885488

beaustates <- c("State 0","State 1", "State 2", "State 3", "State 4")

beaucoeff <- c(0, -0.3979885, -0.5518057, -0.8339985, -1.3146111)

esw <- vector(mode="numeric", length=length(beaustates))

for (i in 1:length(beaustates)) {

  esw[i] <- integrate(gz, 0, maxz,

                      beta=betaval, sigintercept = sigint, sigcoef = beaucoeff[i],  DistWin=FALSE)$value

  if(i==1) {

    plot(seq(0,maxz),

         gz(beta=betaval, sigintercept = sigint, sigcoef = beaucoeff[i], z=seq(0,maxz), DistWin=FALSE),

         lwd=3, type="l", col=i, main="Detection functions by sea state",

         xlab="Distance", ylab="Detection probability", ylim=c(0,1))

  } else {

    lines(seq(0,maxz),

          gz(beta=betaval, sigintercept = sigint, sigcoef = beaucoeff[i], z=seq(0,maxz), DistWin=FALSE),

          lwd=3, type="l", col=i)

  }

}

msg <- paste(beaustates, title="Estimated ESW", round(esw, 1))

legend("topright", legend=msg, lwd=2, col=1:5)

---

From: distance...@googlegroups.com <distance...@googlegroups.com> on behalf of Tim Awbery <tim.a...@sams.ac.uk>

Sent: 12 May 2023 16:51

To view this discussion on the web visit https://groups.google.com/d/msgid/distance-sampling/e4413ffc-a4c4-4ede-aa24-a6bfaa3c1cb7n%40googlegroups.com.

[Tim Awbery]

Hi Eric,

Thank you very much for getting back to me and thank you to yourself and the team for  working so hard to provide such a valuable resource.

Best,

Tim

[Tiffany Goh]

Hi Eric,

Just piggy-backing on this post!

I understand that the above codes extracts the esw for each factor of the covariate only when using hazard rate (thanks for the codes btw, very useful). How would you tweak the code if it was half normal?

Many thanks,

Tiffany

[Eric Rexstad]

Tiffany

I've extended the previous function to incorporate both the half normal and hazard rate key functions.  Attached

---

From: distance...@googlegroups.com <distance...@googlegroups.com> on behalf of Tiffany Goh <tiffan...@gmail.com>
Sent: 19 January 2024 18:43

To view this discussion on the web visit https://groups.google.com/d/msgid/distance-sampling/b516ef20-036d-4505-8ade-36588689401fn%40googlegroups.com.