Skip to content

Slope credible intervals do not cover slope estimate #45

@stweb75

Description

@stweb75

I have found the slope sign probabilities (slpSgnPosPr ...) easier to interpret than the slope Credible Interval (slpCI).
Nonetheless, the slpCI is often wrong in that the slope estimate (slp) is outside of the CI.

## Example 1; trend only model
library(Rbeast)
data(Nile)
out <- beast(Nile, season="none", dump.ci=TRUE)
plot(out)
# trend CI looks OK
plot(out[[7]][13]$Y, type="l", ylim=c(700,1200))
  lines(out[[7]][15]$CI[,1], lty=3)
  lines(out[[7]][15]$CI[,2], lty=3)
# slope outside slope CI at a few indices 
plot(out[[7]][16]$slp, type="l", ylim=c(-10,10))
  lines(out[[7]][18]$slpCI[,1], lty=3)
  lines(out[[7]][18]$slpCI[,2], lty=3)
which(out[[7]][18]$slpCI[,2] > out[[7]][16]$slp)  
which(out[[7]][18]$slpCI[,1] < out[[7]][16]$slp)'
## Example 2; trend + season model
data(googletrend_beach)
out <- beast(googletrend_beach, dump.ci=TRUE)
plot(out)
# trend CI looks OK
plot(out[[7]][13]$Y, type="l", ylim=c(55,75))
  lines(out[[7]][15]$CI[,1], lty=3)
  lines(out[[7]][15]$CI[,2], lty=3)
# season CI looks OK
plot(out[[8]][13]$Y, type="l", ylim=c(-20,30))
  lines(out[[8]][15]$CI[,1], lty=3)
  lines(out[[8]][15]$CI[,2], lty=3)
# season Y within CI  
which(out[[8]][15]$CI[,2] > out[[8]][13]$Y)  
which(out[[8]][15]$CI[,1] < out[[8]][13]$Y) 
# slope frequently outside slope CI
plot(out[[7]][16]$slp, type="l", ylim=c(-10,10))
  lines(out[[7]][18]$slpCI[,1], lty=3)
  lines(out[[7]][18]$slpCI[,2], lty=3)
which(out[[7]][18]$slpCI[,2] > out[[7]][16]$slp)  
which(out[[7]][18]$slpCI[,1] < out[[7]][16]$slp)

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions