4
$\begingroup$

EDIT: Thanks for the help so far. I have updated the question based on further work.

I am interested in finding any difference in response of counts of two competing species to dryness and elevation. Both are caught in the same type of trap, but trap placement (bushes vs buildings) probably biases toward one species or the other, so each trap placement is rated according to indoors/outdoors on a scale of 1-3. A trap does not move once it is placed. Traps are nested in sites, dryness is measured at each site/visit combination, and elevation is of course constant at each site.

My data includes the number of each species caught in 4-12 traps, at 8-12 sites during 6-8 visits to each site. The visits and sites are far enough spaced that I assume they are independent. Most data are zeros, especially for Sp_A which is rare. Below is the structure of the data.

R> str(mydata)
'data.frame':   300 obs. of  8 variables:
     $ count    : num  0 5 1 0 1 1 0 0 0 0 ...
     $ species  : Factor w/ 2 levels "a","b": 1 1 1 1 1 1 1 1 1 1 ...
     $ elevation: int  1 1 1 1 1 1 1 1 1 1 ...
     $ dryness  : num  0.179 0.179 0.179 0.179 0.179 ...
     $ site     : Factor w/ 5 levels "a","b","c","d",..: 1 1 1 1 1 1 1 1 1 1 ...
     $ trap     : Factor w/ 50 levels "MT10a","MT10b",..: 6 11 16 21 26 31 36 41 46 1 ...
     $ visit    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
     $ in_out   : num  1 1 1 1 1 1 1 1 1 1 ...

As an initial attempt, I had no trouble using R function lmer to model Sp_A presence/absence with Sp_B as a covariate: (SpA>0) ~ Sp_B + elevation + (1|site/trap) + (1|visit)

However on suggestion from R-sig-mixed-models list, and also "ils" from this site, I have reshaped the data as seen above, and tried to model count with species as a factor in the model. I have used two different packages (lme4 and glmmADMB). The functions handle the data ok in simple models, until I add the factor "species." Below are the R code and error messages. These same errors also happen when I use simpler dummy data. Any ideas on syntax or other packages?

require(lme4)
glmera1 <- glmer(count~elevation*species*in_out 
                 + (1|visit) + (1|site/trap), 
                 data=mydata, family="poisson")
Error in asMethod(object) : matrix is not symmetric [1,2]
In addition: Warning messages:
1: In mer_finalize(ans) :
  Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 432
2: In mer_finalize(ans) :
  Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 432
3: In mer_finalize(ans) : singular convergence (7)

require(glmmADMB)
admba1 <- glmm.admb(count~elevation*species*in_out, random=~visit, 
                    group="site", data=mydata, family="nbinom")
Error in glmm.admb(count ~ elevation * species * in_out, random = ~visit,  : 
  The function maximizer failed
In addition: Warning message:
running command './nbmm -maxfn 500 ' had status 1
Improve this question
$\endgroup$
7
  • $\begingroup$ I'm not aware of any question that covers "exactly the same ground." While these criteria are sometimes negative, my question is, have I done something wrong? $\endgroup$
    – J. Win.
    Apr 7, 2011 at 22:52
  • 2
    $\begingroup$ @J. Winchester I'm afraid your question still appears to focus on AIC, and responses to the related question on AIC apply as well here. However, I feel you are also wondering whether your lmer syntax allows you to test the hypothesis formulated in your first paragraph. If this is the case, you can simply edit the title and your question and focus on this issue, rather than AIC. You might also consider adding the output of summary(lmera.1). Finally, the outcome you consider is a ratio of counts (% Sp_A): Is it usual in your field to treat this as a continuous response variable? $\endgroup$
    – chl
    Apr 8, 2011 at 19:30
  • $\begingroup$ @chl. Thank you. I have rewritten the question and title; regarding your last sentence, I cannot say that it is usual, but I do not have a better idea. $\endgroup$
    – J. Win.
    Apr 12, 2011 at 8:46
  • 1
    $\begingroup$ @J. Winchester: It might be a crazy idea, but why not treat abundace as the dependent variable and add species_identity as a factor into the model? Then, when you have estimated the magnitude of species_identity you will be able to compare the abundance of both species in relative terms. This way you won't need the model comparison stage and things will be much more simple. I also don't see the purpose of visit_day as an estimable parameter; the species abundance will clearly depend on it as an exposure parameter, but then how one will interpret the dryness * visit_day interaction? $\endgroup$
    – ils
    Apr 12, 2011 at 11:34
  • $\begingroup$ @ils Thank you, I have tried that and updated the question. I wonder what is the difference between your suggestion, and using lmer(cbind(sp_a, sp_b)~elevation*dryness ... ) which is a format used at the bottom of the lmer help page. $\endgroup$
    – J. Win.
    Apr 20, 2011 at 23:38

2 Answers 2

Reset to default
3
$\begingroup$

Check out this question: Negative values for AICc (corrected Akaike Information Criterion)

The same explanation holds for BIC.

Improve this answer
$\endgroup$
1
$\begingroup$

Try the following packages on cran: Rcapture and mra. The easiest method to use in R for a log linear model, open closed or robust, M0, Mh, Mt, or Mth is well detailed here: http://www.jstor.org/pss/3695557 and that is what Rcapture is based on. You forgot to specify in your lmer call that it should be a count response (family = poisson) and your Y should be the abundance associated with a particular capture history. Extracting the estimated abundances, survival probilities, etc, is well detailed in the jstor link.

However, the premier tool people use for this is MARK (http://www.phidot.org/software/mark/docs/book/), and there is a bridge package not on cran called RMark (http://www.phidot.org/software/mark/rmark/RMarkWorkshopNotes.pdf). RMark is the premier tool for this sort of thing and the one which certainly has the largest community. It has its own support forum at phidot.org.

Also, google around for Dr. Schofield's BUGS code for this problem; it is online somewhere, I think in the dissertation.

Btw, if you are interested in a job in population ecology and R programming, we're hiring: https://performancemanager4.successfactors.com/career?company=comScore&career_job_req_id=721&career_ns=job_listing&navBarLevel=JOB_SEARCH

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.