EliGurarie / cyclomort

R package for survival modeling with periodic mortality function
2 stars 0 forks source link

Axes inconsistent when plotting cmfactorfit #49

Open coverton-usgs opened 3 years ago

coverton-usgs commented 3 years ago

Wonderful package! I am so glad to have just found it.

I am working on contrasting seasonal survival patterns of two populations of birds. When summarizing the cmfactorfit object there is no significant evidence of a difference between groups which I want to show with a graphic:

Summary table comparing factorial seasonal survival model with 1 seasons.

Formula: event ~ classfactor
[[1]]
[[1]]$test
              logLike k    AIC  LRT p.value
null        -43.12865 3 92.257 0.22   0.974
classfactor -43.01842 6 98.037             

[[1]]$null
[[1]]$null$model
[1] "Multi-seasonal hazard function fit with 1 seasons with periodicity 365.\n\n"

[[1]]$null$estimates
[[1]]$null$estimates$meanhazard
     estimate      CI.low     CI.high           se
1 0.004374442 0.003426811 0.005322072 0.0004738155

[[1]]$null$estimates$`season 1`
  parameter season estimate    CI.low  CI.high
1      peak      1 299.1501 263.54784 334.7523
2  duration      1 133.8589  89.98236 161.3743

[[1]]$null$analysis
[1] "Log-likelihood: -43.1287; AIC: 92.2573"

[[2]]
[[2]]$R
[[2]]$R$model
[1] "Multi-seasonal hazard function fit with 1 seasons with periodicity 365.\n\n"

[[2]]$R$estimates
[[2]]$R$estimates$meanhazard
     estimate      CI.low     CI.high          se
1 0.004353411 0.003389777 0.005317045 0.000481817

[[2]]$R$estimates$`season 1`
  parameter season estimate    CI.low  CI.high
1      peak      1 298.6343 262.00630 335.2623
2  duration      1 133.6692  89.89374 161.2140

[[2]]$R$analysis
[1] "Log-likelihood: -41.6957; AIC: 89.3914"

[[2]]$T
[[2]]$T$model
[1] "Multi-seasonal hazard function fit with 1 seasons with periodicity 365.\n\n"

[[2]]$T$estimates
[[2]]$T$estimates$meanhazard
     estimate        CI.low     CI.high          se
1 0.004508746 -0.0007287932 0.009746284 0.002618769

[[2]]$T$estimates$`season 1`
  parameter season  estimate       CI.low  CI.high
1      peak      1  91.95443 -130.0376013 313.9465
2  duration      1 133.75070    0.5225195 182.4227

[[2]]$T$analysis
[1] "Log-likelihood: -1.3227; AIC: 8.6454"

However, the resulting plot places the "Null" model centered on a y-axis scaled from 0-1 using fit = "both" or fit = "null" The "Alt" plots appear to be placed approrpiately for the 0-1 y axis when plotting fit = "both"

Rplot02

However when plotting fit = "alt" the y axis is only scaled for first factor.
Rplot01

The above graphic should have lines that almost completely overlap (with factor #2's confidence intervals covering the whole range of the y-axis

If I plot with a specified axis range it also doesn't match the figure below is with:

plot(cyclo_fit_gt15d_fact, fit = "alt", ymax = 0.1)
title(main=paste0("sequestered days = ", sequestereddays))

Rplot

Sorry I know its bad form to post two issues at once but I didn't realize the null model wasn't scaled with appropriate ylims labels until writing this.

EliGurarie commented 3 years ago

Glad the package is (potentially) useful for your birds!

It's a bit tricky to tell what's going on with the axes without a minimal working example - and I have very little bandwidth at the moment to play with this too much. But if you just want to illustrate that the two fits are very similar (as the stats certainly show) you can just fit the separate models and use the add = TRUE tag, e.g.

fit1 <- fit_cyclomort(T.morts.pop1, n.seasons = 1)
fit2 <- fit_cyclomort(T.morts.pop2, n.seasons = 1)

plot(fit1, hazcolor = "blue")
plot(fit2, add = TRUE, hazcolor = "red")

That's pretty simple and should be robust.

Also - if you want to tweak the plots in any other way, you can check out the guts of the plot.cmfit and the - much fussier - plot.cmfactorfit functions.

Maybe if you mess around, you can see where there might be an error?