model-averaging with multiple covariates

[PaulS]

Hi,
I have multi-season data that I am analyzing with colext models. I haven't used model-averaging much in the past and I'm really struggling with it, and would appreciate some guidance.  Say for example I have 8 models (listed below) that each have different combinations of covariates (2 continuous, 2 categorical) influencing psi.  For simplicity, I've left the col, ext, and p parameters as constant here.

#Sites: 31
#Primary sampling periods: 5
#Secondary sampling occasions: 3

#Covariates influencing psi:
#elev.site=700:1600
#rip.site=0:1.5
#res.site='low', 'none'
#burn.site='C', 'E', 'I'

m1 <- colext(~burn.site, ~1, ~1, ~1, umf)
m2 <- colext(~res.site, ~1, ~1, ~1, umf)
m3 <- colext(~burn.site + res.site, ~1, ~1, ~1, umf)
m4 <- colext(~1, ~1, ~1, ~1, umf)
m5 <- colext(~burn.site + scale(elev.site), ~1, ~1, ~1, umf)
m6 <- colext(~res.site + scale(elev.site), ~1, ~1, ~1, umf)
m7 <- colext(~burn.site + scale(rip.site), ~1, ~1, ~1, umf)
m8 <- colext(~res.site + scale(rip.site), ~1, ~1, ~1, umf)

#I then create a fitList and run the model selection:

fmList <- fitList(m1=m1, m2=m2, m3=m3, m4=m4, m5=m5, m6=m6, m7=m7, m8=m8)
modSel(fmList)

#Once I have the model selection output, I am wondering how you get model-average predictions for each covariate. Because I have 4 different covariates in my set of candidate models, I need to include all 4 of them in the dataframe, correct?

data.predict <- data.frame(elev.site=seq(700, 1600, length=31), rip.site=seq(0, 1.5, length=31),
                           burn.site=rep(c("C", "E", "I"), length=31),
                           res.site=rep(c("low", "none"), length=31))
                          
mod.pred <- predict(fmList, type='psi', newdata=data.predict, appendData=TRUE)
mod.pred

#However, this gives all possible combinations and doesn't really allow you to tease apart each covariate.  From reading some of the other posts and going through the vignette (2012), it seems like I should set the continuous variables to 0 and/or set one of the categorical variables to a particular level.  So, I tried something like this:

data.predict <- data.frame(elev.site=0, rip.site=0,
                           burn.site=rep(c("C", "E", "I")),
                           res.site=factor('low', levels=c("low", "none")))
mod.pred <- predict(fmList, type='psi', newdata=data.predict, appendData=TRUE)
mod.pred

#But, I get this error:

Error in X[i, ] : subscript out of bounds

Any insights on how to get these estimates so I can present them in a figure, or should I compare a simpler set of models?  I would really appreciate your help.

Thanks.

[Richard Chandler]

Hi Paul,

I don't know what is causing that error. Can you save() the object "fmList" and send it to me? So that the file isn't too big, you could do this:

save(fmList, file="fmList.gzip")

and email me the file "fmList.gzip"

Richard

--
You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[PaulS]

Hi Richard,

Thanks for your response. I've been fiddling around with this today and changed the code to basically give it 0 values, but in a sequence:

data.predict <- data.frame(elev.site=seq(0, 0.00001, length=3), rip.site=seq(0, 0.00001, length=3)
                           res.site=factor("none", levels=c("low", "none")),
                           burn.site=rep(c('C', 'E', 'I'), length=3))

mod.pred <- predict(fmList, type='psi', newdata=data.predict, appendData=TRUE)
mod.pred

And this worked.  Do you see any issues with simply proceeding this way if it seems to be working? 

If not, I can work on getting you the file you requested.

Thanks.

[Richard Chandler]

Oh right, you can't do: scale(0). If you standardize your covariates before putting them in the unmarkedFrame, this problem will go away.

Richard

[Whitney Wiest]

Hi Paul and Richard,

I'm getting the same error message (Error in X[i, ] : subscript out of bounds) when I try to predict detection probability for all levels of one covariate, while holding the other model covariates constant.  What I'm really trying to do eventually is predict detection for each level of each factor when you have missing data for the other factors.  Is this possible?  Any insight is greatly appreciated!  Thanks! Whitney


##Some code

summary(will.occu)
unmarkedFrame Object

59 sites
Maximum number of observations per site: 8
Mean number of observations per site: 8
Sites with at least one detection: 51

Tabulation of y observations:
   0    1    2    3    4    5    6    7   12 <NA>
 276   80   53   24   21   11    5    1    1    0

Observation-level covariates:
 noise   tide    sky     wind    date    observer
 0:142   1: 20   0:217   0: 42   1:113   1: 44  
 1:247   2:106   1:115   1:142   2:294   2: 48  
 2: 80   3: 82   2:122   2:178   3: 65   3:104  
 3:  3   4: 27   4:  9   3: 91           4: 44  
         5:119   5:  8   4: 19           5: 40  
         6:118   8:  1                   6: 71  
                                         7:121  
> str(will.occu)
Formal class 'unmarkedFrameOccu' [package "unmarked"] with 5 slots
  ..@ y       : num [1:59, 1:8] 0 0 0 0 0 0 4 1 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:59] "178" "179" "180" "181" ...
  .. .. ..$ : chr [1:8] "count1" "count2" "count3" "count4" ...
  ..@ obsCovs :'data.frame':    472 obs. of  6 variables:
  .. ..$ noise   : Factor w/ 4 levels "0","1","2","3": 3 3 3 3 3 2 2 2 3 3 ...
  .. ..$ tide    : Factor w/ 6 levels "1","2","3","4",..: 2 5 3 6 6 6 6 5 6 5 ...
  .. ..$ sky     : Factor w/ 6 levels "0","1","2","4",..: 1 3 1 3 1 3 1 3 1 3 ...
  .. ..$ wind    : Factor w/ 5 levels "0","1","2","3",..: 4 4 3 3 4 2 2 2 3 3 ...
  .. ..$ date    : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 2 1 1 ...
  .. ..$ observer: Factor w/ 7 levels "1","2","3","4",..: 5 5 5 5 5 5 5 5 5 5 ...
  ..@ siteCovs: NULL
  ..@ mapInfo : NULL
  ..@ obsToY  : num [1:8, 1:8] 1 0 0 0 0 0 0 0 0 1 ...
> will.occu.8global <- occu(~noise+sky+tide+wind+date+observer ~1, will.occu)
Warning message:
In truncateToBinary(y) :
  Some observations were > 1.  These were truncated to 1.
> will.occu.8global

Call:
occu(formula = ~noise + sky + tide + wind + date + observer ~
    1, data = will.occu)

Occupancy:
 Estimate    SE    z  P(>|z|)
     1.93 0.406 4.75 2.07e-06

Detection:
            Estimate      SE       z  P(>|z|)
(Intercept)   1.6736   0.875  1.9124 5.58e-02
noise1       -0.3323   0.322 -1.0316 3.02e-01
noise2       -0.3098   0.426 -0.7280 4.67e-01
noise3       -2.4282   1.582 -1.5351 1.25e-01
sky1          0.2527   0.308  0.8199 4.12e-01
sky2          0.8638   0.301  2.8705 4.10e-03
sky4         -0.3687   0.906 -0.4068 6.84e-01
sky5          1.3024   1.083  1.2028 2.29e-01
sky8         12.9472 504.935  0.0256 9.80e-01
tide2        -0.0218   0.595 -0.0366 9.71e-01
tide3        -0.6554   0.609 -1.0760 2.82e-01
tide4        -0.6808   0.729 -0.9334 3.51e-01
tide5        -0.1950   0.594 -0.3284 7.43e-01
tide6        -0.6263   0.593 -1.0563 2.91e-01
wind1         0.2016   0.470  0.4292 6.68e-01
wind2         0.2194   0.469  0.4680 6.40e-01
wind3        -0.0779   0.521 -0.1496 8.81e-01
wind4         0.9379   0.731  1.2830 1.99e-01
date2        -0.2036   0.277 -0.7347 4.63e-01
date3        -2.4360   0.496 -4.9125 8.99e-07
observer2    -1.2100   0.565 -2.1411 3.23e-02
observer3    -1.2184   0.537 -2.2684 2.33e-02
observer4    -1.7707   0.627 -2.8241 4.74e-03
observer5    -0.3389   0.591 -0.5735 5.66e-01
observer6    -1.8285   0.519 -3.5213 4.29e-04
observer7    -1.4173   0.491 -2.8840 3.93e-03

AIC: 577.3579

##Is there a way to set sky - observer as "missing values" and calculate detection for all levels of noise?

> will8.noise <- data.frame(noise=factor(c("0","1","2","3"),levels=c("0","1","2","3")), sky=factor("0",levels=c("0","1","2","4","5","8")), tide=factor("1",levels=c("1","2","3","4","5","6")), wind=factor("0",levels=c("1","2","3","4")), date=factor("1",levels=c("1","2","3")), observer=factor("1",levels=c("1","2","3","4","5","6","7")))

> predict(will.occu.8global, type = "det", newdata = will8.noise, appendData = TRUE)

Error in X[i, ] : subscript out of bounds

[Richard Chandler]

Hi Whitney,

It looks like you left out the "0" level for wind:

wind=factor("0",levels=c("1","2","3","4"))

Richard

Richard Chandler

Assistant Professor

Warnell School of Forestry and Natural Resources

University of Georgia

706-542-5815

chandlerlab.wordpress.com

[Whitney Wiest]

Richard,

I sure did.  Thanks, that fixed the original error message.  Is there a way to code "sky" through "observer" as missing data values and then estimate detection for each level of noise, i.e., to mimic when a surveyor "forgets" to record some Obs cov values?

Thanks,
Whitney

> will8.noise <- data.frame(noise=factor(c("0","1","2","3"),levels=c("0","1","2","3")), sky=factor("0",levels=c("0","1","2","4","5","8")), tide=factor("1",levels=c("1","2","3","4","5","6")), wind=factor("0",levels=c("0","1","2","3","4")), date=factor("1",levels=c("1","2","3")), observer=factor("1",levels=c("1","2","3","4","5","6","7")))

> predict(will.occu.8global, type = "det", newdata = will8.noise, appendData = TRUE)

  Predicted        SE      lower     upper                  noise sky tide wind date observer
1 0.8420611 0.1163933 0.48958421 0.9673577     0   0    1    0    1        1
2 0.7927064 0.1352645 0.43240327 0.9504845     1   0    1    0    1        1
3 0.7963816 0.1396321 0.41973327 0.9548485     2   0    1    0    1        1
4 0.3198314 0.3571927 0.01847395 0.9215542     3   0    1    0    1        1

[Richard Chandler]

There is no easy way to do what you want to do. I think you would need to choose priors for your two covariates and then average over them, using either the bootstrap or Bayesian methods. 

[Whitney Wiest]

Ok, thanks for the head's up!

Whitney

Whitney A. Wiest
250 Townsend Hall
Department of Entomology and Wildlife Ecology
University of Delaware
Newark, DE 19717
wwi...@udel.edu

U.S. Fish and Wildlife Service
Indefinite Pathways Intern, Biological Sciences
 

P Please consider the environment before you print this email

--
You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/f4IroYMN1Jg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+u...@googlegroups.com.