Ah yes. Forgot about relevel().  That would  be simpler.

Cheers,
Bert
Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Thu, Feb 25, 2016 at 1:19 PM, Marc Schwartz <marc_schwa...@me.com> wrote:
> Hi,
>
> Well...as I understand the question:
>
> If Katharine wants to use treatment contrasts, which are the default, the 
> easiest thing to do may be use to ?relevel to specify the reference level for 
> the IV factor.
>
> However, as Bert notes, there are other contrasts that can be used, which 
> would affect the interpretation of the comparisons of the factor levels and 
> the help pages that he references would cover a high level review of the 
> options, with more detail provided in the reference therein.
>
> Also, if there is an excess of zeroes, you may need to consider the use of a 
> zero inflated model and the 'pscl' package on CRAN is worthy of 
> consideration, along with its vignette on count regression models:
>
> https://cran.r-project.org/web/packages/pscl/
> https://cran.r-project.org/web/packages/pscl/vignettes/countreg.pdf
>
> The ?vuong test in that package can also be helpful for comparing 
> zero-inflated models with non zero-inflated models.
>
> Regards,
>
> Marc Schwartz
>
>
>> On Feb 25, 2016, at 2:48 PM, Bert Gunter <bgunter.4...@gmail.com> wrote:
>>
>> You can re-set the contrasts for the factor, though whether this is
>> "easier" is a matter of personal preference.
>>
>> See ?C or ?contrasts, which I understand to be alternative ways of
>> doing the same thing (and would appreciate correction is this is
>> wrong). Or this can be done through the "contrasts" argument of
>> glm.nb.
>>
>> Cheers,
>> Bert
>>
>>
>>
>>
>> Bert Gunter
>>
>> "The trouble with having an open mind is that people keep coming along
>> and sticking things into it."
>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>>
>>
>> On Thu, Feb 25, 2016 at 10:53 AM, Katharine Miller - NOAA Federal
>> <katharine.mil...@noaa.gov> wrote:
>>> Hi,
>>> I am using the glm.nb function to evaluate differences in the catch of
>>> individual species of fish across three river tributaries.  The dependent
>>> variable is catch per unit effort (CPUE), and the independent variable is
>>> the the tributary (Trib_cat).  CPUE is derived from the fish counts divided
>>> by the effort, so the response is not a count per se, but I think the
>>> negative binomial is appropriate because of the large numbers of zeros in
>>> the dependent variable. Trib_cat is a column with 1, 2, or 3 depending on
>>> which tributary it is representing.
>>>
>>> I have the following code:
>>>
>>> dataS<-read.table("OtherSpecies.txt", sep="", header=T)
>>>
>>> dataS <- within(dataS, {
>>>  Trib_cat <- factor(Trib_cat, levels = 1:3, labels = c("MM", "NM", "SM"))
>>>  Year<-factor(Year)
>>> })
>>>
>>> ## separate out the species of interest
>>> coregonid<-subset(dataS, Spec_age=="Coregonid")
>>>
>>> fit.CT<-glm.nb(CPUE ~ Trib_cat, data=coregonid, link = log)
>>> summary(fit.CT)
>>>
>>> The result comes out like this:
>>> Call:
>>> glm.nb(formula = CPUE ~ Trib_cat, data = coregonid4, init.theta =
>>> 0.1723775759,
>>>    link = log)
>>>
>>> Deviance Residuals:
>>>    Min       1Q   Median       3Q      Max
>>> -0.9239  -0.8055  -0.6687  -0.6286   2.9020
>>>
>>> Coefficients:
>>>            Estimate Std. Error z value Pr(>|z|)
>>> (Intercept)  -0.7805     0.2497  -3.125  0.00178 **
>>> Trib_catNM   -0.2140     0.3485  -0.614  0.53921
>>> Trib_catSM    0.7393     0.3296   2.243  0.02488 *
>>> ---
>>>
>>> This gives me the difference in CPUE between the NM and SM tributaries
>>> compared to the MM tributary  (acting here as the reference group).  What I
>>> need is to compare all of the tributaries with all the others - so I need
>>> to run the model twice with a different reference group on the second run.
>>>  I can do this by changing the numbering of the Trib_cat field  in the
>>> underlying database (e.g. changing MM from 1 to 3, SM from 2 to 1, etc) and
>>> re-running the model, but I, as I have a number of species to do this for,
>>> I was wondering if there was an easier way to specify which group to use as
>>> the reference group when calling the model.
>>>
>>> Thanks for any help.
>

______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to