1
$\begingroup$

I've got what I think is a fairly basic problem. I'm not sure if it is a conceptual question or software question, but I'm fairly new to using R and these kinds of stats, so it could be either or both.

I'm working on a practice problem with the occurence of a certain bird in fields with various combinations of two vegetation types (hydric or mesic) and two perscribed burn treatments (burned or not burned). The data set oringally contained counts of birds in each field, but since detections of the birds are fairly rare (i.e., there are a lot of 0's in the counts), we are using presence/absence of the bird as our response and using a generalized linear model with a binomial distribution. The data (which includes columns of data for other problems) look something like:

> birddata
   vegtype d_veg       burn d_burn birdsum birdpa  offset     loff    lnoff  forb rcdom
1   hydric     0 not burned      1       0      0 25.1328 1.400241 3.224174  0.10  5.20
2   hydric     0     burned      0       7      1 25.1328 1.400241 3.224174  6.55  5.20
3   hydric     0     burned      0       0      0 25.1328 1.400241 3.224174  6.40  4.60
4   hydric     0     burned      0       3      1 25.1328 1.400241 3.224174 13.35  4.45
5   hydric     0 not burned      1       0      0 25.1328 1.400241 3.224174 11.70  4.20
6   hydric     0     burned      0       0      0 25.1328 1.400241 3.224174 19.10  0.80
7   hydric     0     burned      0       0      0 25.1328 1.400241 3.224174  1.90  5.10
8    mesic     1     burned      0       6      1 25.1328 1.400241 3.224174 12.95  0.05
...
49   mesic     1 not burned      1       2      1 25.1328 1.400241 3.224174 24.40  0.05
50   mesic     1 not burned      1      14      1 25.1328 1.400241 3.224174  4.10  1.05

and the model we're using is:

> glm.birdpa2 <- glm(birdpa ~ burn + vegtype, family = "binomial", data = birddata)

What I've been trying to figure out is how to use R to get the probability of detecting a bird in each of the vegetation types while holding the burn treatment at their mean value. I think what this problem is trying to get at is: how could you tell someone what the overall chance of seeing one of these birds is in a field of either particular vegetation type? I think in SAS I would be using an estimate statement to do ths, but I don't know how to do this in R. I've been trying to use the predict and fitted functions but everything I've tried so far has seemed to be a dead end. Does anyone have some tips for how I'd go about doing this?

Improve this question
$\endgroup$
2
  • $\begingroup$ When you say "holding the burn treatment at their mean value" do you mean "at the overall mean across both vegetation types" or "at the mean value within vegetation type"? $\endgroup$
    – Glen_b
    Commented Oct 18, 2014 at 22:18
  • $\begingroup$ Glen_b, that's a good question, and the problem is not very clear (the bold text above is almost verbatim). I think "at the mean value within vegetation types" would be more appropriate though. $\endgroup$
    – chrysoparia
    Commented Oct 19, 2014 at 14:43

1 Answer 1

Reset to default
1
$\begingroup$

I would use the lsmeans package to do these predicted probabilities. Load the package with library(lsmeans), and then execute the following to get predicted probabilities of detecting a bird at each vegetation type, averaged over the levels of burn treatment:

lsmeans(glm.birdpa2,
        ~ vegtype,
        type="response",
        data=birddata)

The following will give you predicted probabilities of detecting a bird at each vegetation type within each level of burn treatment:

lsmeans(glm.birdpa2,
        ~ vegtype | burn,
        type="response",
        data=birddata)

The type="response" argument in the above tells lsmeans to back-transform the results, so the predictions will be probabilities. Leaving that argument out, lsmeans will provide log-odds rather than probabilities.

Improve this answer
$\endgroup$

Your Answer

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

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