Counterintuitive/unintended behavior of get.real?

[Daniel Barton]

Hi All (I guess mainly Chris, Andy, Dan),

I'm finally working on publishing this small mammal thing I talked to Andy about using oSCR for a couple years ago.  I'm finding what I thought was an odd issue, but maybe it's user error. 

If I:

m7  <- oSCR.fit(model=list(D~cov+session,p0~b,sig~sex), scrFrame=lassicFrame, ssDF=lassicSS, trimS=150, multicatch = TRUE, sexmod="session")

then:

get.real(model=m7, type="dens")

(just an example...) then this doesn't yield sex-specific density estimates - no matter what newdata I do or don't pass to it.  I see there is a seemingly undocumented argument to function get.real in the source, N_sex=T/F, but this appears to have no effect.  I can produce sex-specific density estimates "by hand" but that's a little annoying and it doesn't seem like this is the intended behavior of get.real() when using it on an oSCR.fit model object.  Am I just doing something hare-brained, or should I dig into the source?

I'm digging the oSCR package still, so thanks for all this.

THANKS,

Dan

--

--
Dan Barton

Note: this is not my official Humboldt State email.

[Chris Sutherland]

Hello Dan,

Couple of questions, what is the output when you type:

print.oSCR.fit(m7)

Second, what happens when you do:

ns <- length(m7$ssDF)

nd <- data.frame(cov = rep(seq(min(cov), max(cov), length=30),ns), session=factor(rep(1:ns, each=30)))

get.real(model=m7, newdata=nd, type="dens")

(On my phone here, so code may have wee bugs, but hopefully you get what I mean)

Chris

--
You received this message because you are subscribed to the Google Groups "oSCR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to oscr_package...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/oscr_package/CAGSCbgpR1GZJenxTY31t-VCgKbNiN6d0a0tvenWxh4f03GDr-Q%40mail.gmail.com.

[Daniel Barton]

Thanks Chris!  Here's output for each of those:

print.oSCR.fit(m7)  

Model:  D ~ cov + session p0 ~ b sig ~ sex
 Run time:  1139.412  minutes
 AIC:  7160.198
 
Summary table:
                     Estimate    SE       z P(>|z|)
p0.(Intercept)         -2.850 0.091 -31.445   0.000
p.behav                 1.851 0.159  11.607   0.000
sig.(Intercept)         3.518 0.037  93.963   0.000
sig.sexmale            -0.613 0.080  -7.652   0.000
d0.(Intercept)         -2.774 0.263 -10.552   0.000
d.beta.covburnfor      -0.615 0.190  -3.234   0.001
d.beta.covchaparral     1.032 0.191   5.406   0.000
d.beta.covforest        0.286 0.318   0.900   0.368
d.beta.covserpentine   -0.939 0.373  -2.519   0.012
d.beta.session2         0.397 0.198   2.007   0.045
d.beta.session3         1.060 0.229   4.617   0.000
d.beta.session4         1.041 0.232   4.482   0.000
d.beta.session5         0.755 0.246   3.072   0.002
d.beta.session6         0.822 0.240   3.424   0.001
d.beta.session7         0.832 0.239   3.476   0.001
psi1                   -0.465 0.395  -1.176   0.240
psi2                    0.100 0.244   0.408   0.683
psi3                   -1.024 0.299  -3.421   0.001
psi4                   -2.596 0.598  -4.340   0.000
psi5                   -0.668 0.354  -1.888   0.059
psi6                   -1.376 0.415  -3.315   0.001
psi7                   -1.577 0.442  -3.569   0.000
*Density intercept is log(individuals per pixel)
  Nhat(state-space) = exp(d0.)*nrow(ssDF)
  (caution is warranted when model contains density covariates)

[note that "cov" is actually a factor because my naming conventions are trash]

ns <- length(m7$ssDF)
nd <- data.frame(cov = rep(c("burnchap","burnfor","chaparral","forest","serpentine"),ns), session=factor(rep(1:ns, each=5)))

get.real(model=m7, newdata=nd, type="dens")

yields

     estimate          se        lwr        upr        cov session
1  0.06240166 0.016406262 0.03727343 0.10447032   burnchap       1
2  0.03374763 0.009144829 0.01984189 0.05739889    burnfor       1
3  0.17519700 0.033282007 0.12073182 0.25423280  chaparral       1
4  0.08310164 0.020648084 0.05106337 0.13524141     forest       1
5  0.02439653 0.010549099 0.01045347 0.05693715 serpentine       1
6  0.09278915 0.021519425 0.05889590 0.14618720   burnchap       2
7  0.05018158 0.012119383 0.03125837 0.08056053    burnfor       2
8  0.26051199 0.037203379 0.19690942 0.34465846  chaparral       2
9  0.12356932 0.027581477 0.07978356 0.19138499     forest       2
10 0.03627681 0.015073632 0.01606697 0.08190762 serpentine       2
11 0.18003426 0.027527952 0.13341378 0.24294594   burnchap       3
12 0.09736487 0.017095791 0.06901468 0.13736089    burnfor       3
13 0.50545871 0.079139188 0.37188736 0.68700508  chaparral       3
14 0.23975552 0.071144388 0.13402344 0.42890039     forest       3
15 0.07038613 0.025915909 0.03420362 0.14484452 serpentine       3
16 0.17680024 0.027206251 0.13076648 0.23903927   burnchap       4
17 0.09561587 0.016757903 0.06781773 0.13480831    burnfor       4
18 0.49637895 0.077827339 0.36504875 0.67495659  chaparral       4
19 0.23544869 0.069812478 0.13167501 0.42100688     forest       4
20 0.06912175 0.025696369 0.03335572 0.14323829 serpentine       4
21 0.13275386 0.024354032 0.09265934 0.19019763   burnchap       5
22 0.07179501 0.014545225 0.04826618 0.10679369    burnfor       5
23 0.37271569 0.068388552 0.26012988 0.53402934  chaparral       5
24 0.17679118 0.055105985 0.09596997 0.32567606     forest       5
25 0.05190140 0.019627075 0.02473330 0.10891208 serpentine       5
26 0.14197827 0.024224709 0.10162135 0.19836215   burnchap       6
27 0.07678369 0.014633730 0.05284949 0.11155710    burnfor       6
28 0.39861386 0.068555157 0.28454880 0.55840337  chaparral       6
29 0.18907553 0.057656838 0.10400746 0.34372107     forest       6
30 0.05550777 0.020850333 0.02658348 0.11590328 serpentine       6
31 0.14336947 0.024118989 0.10309955 0.19936853   burnchap       7
32 0.07753607 0.014603052 0.05360269 0.11215562    burnfor       7
33 0.40251976 0.068391069 0.28850885 0.56158470  chaparral       7
34 0.19092822 0.057993581 0.10527294 0.34627687     forest       7
35 0.05605167 0.021037797 0.02685978 0.11697007 serpentine       7

THANKS.

To view this discussion on the web visit https://groups.google.com/d/msgid/oscr_package/CANjB3vQHvDUfc60_Dtis-r42e5rk2azo7HoA2MXUbVNwy6Sztw%40mail.gmail.com.

[Chris Sutherland]

Hi Dan,

 

Hmm, okay. Any chance you could send me along the m7 object?

 

Chris

To view this discussion on the web visit https://groups.google.com/d/msgid/oscr_package/CAGSCbgo2y6E6yuw4tg-0Gq9sM4R5VczjUnvosBuRNsbp6hDSGQ%40mail.gmail.com.

[Daniel Barton]

Sure, sent off-list... if anyone is following with bated breath I'll follow back up to the list.

Best,

Dan

To view this discussion on the web visit https://groups.google.com/d/msgid/oscr_package/008101d66ea7%243191a270%2494b4e750%24%40gmail.com.

[Chris Sutherland]

Hi Dan and everyone,

 

So this issue was, indeed, related to estimating and predicting session specific sex ratios. A couple of things before the solution:

 

1. There was a naming convention issue where the session specific \psi’s were being named psi1, psi2, …, psiS, whereas we usually use ‘.’s to separate (i.e., psi.1, psi.2, …, psi.S). I have fixed that and pushed the changes to the repo.
2. The main issue here though, was that the prediction code wasn’t complete for predicting session specific-sex ratios using a newdata object (with fix (1), get.real() would use psi.1 for all sessions, which is, of course, wrong). This fix is more involved, and wont be done in the next few weeks. However, there is a workaround which I describe below that amounts to manipulating the ssDF object to look like a newdata object (note that if you don’t provide a newdata object, get.real() makes a prediction for every statespace pixel – which can take a very long time).

 

Okay, the fix (very cool application by the way!):

 

#-----------------------------------------------------------------------------#

# update your oSCR, i pushed a fix to the git repo

library(oSCR)   

library(ggplot2)

 

 

#--------------------------- A HACK ------------------------------------------#

# - because the density covariate is categorical, this works okay

# - note in this specific case, not all categories are present in each

#   session, so I recommend removing those

# - basically we create a new ssDF with one that has a similar structure to

#   the 'newdata' object, but is a list rather than a data frame (the list

#   structure is required, for now, to get session specific psi's)

 

 

levs <- c("burnchap","burnfor","chaparral","forest","serpentine")

ses <- 1:7

m7_forPred <- m7

pred_ss <- m7$ssDF

for(i in ses){

  df <- data.frame(cov = factor(levs, levels = levs),

                   session=factor(i, levels=ses),

                   X=i, Y=i)

  pred_ss[[i]] <- df

}

m7_forPred$ssDF <- pred_ss

 

gr <- do.call(rbind,get.real(model=m7_forPred, type="dens"))

gr$session <- factor(rep(ses,each=10))

gr$cov <- factor(rep(levs,14))

 

 

ggplot(gr,aes(x=session,y=estimate,group=sex, color=sex)) +

  geom_errorbar(aes(ymin=lwr,ymax=upr),position=position_dodge(0.5),width=0) +

  geom_line(linetype=2) +

  geom_point(position = position_dodge(0.5),size=4) +

  facet_wrap(.~cov) +

  theme_bw()

 

 

 

To view this discussion on the web visit https://groups.google.com/d/msgid/oscr_package/CAGSCbgra7%2BzGNENGi%2BtC%2B7OuPhS5NkTNyCzm1kbKJ18g1LSt9A%40mail.gmail.com.

[Daniel Barton]

Thanks a ton, Chris!

To view this discussion on the web visit https://groups.google.com/d/msgid/oscr_package/007701d66f04%241ac70050%24505500f0%24%40gmail.com.