Hello Ben and all,

Thank you Ben, very much for your previous reply.  Yes, stand/plot makes
sense.  After talking with a few statisticians, I got convergence!   The
repeated measure issue is solvable with a dummy variable for each of the 3
years (like a fiscal or school year) - so yr1 = 2010smr+2011wtr, and so
on.  Pellet plots were sampled 6 times, 3 each season, vegetation only
once.  Among year variability can be high.

I'm now attempting to make predictive plots from models.  I used UCLA's
method posted at their IRDA site that has a negative binomial example:

            http://www.ats.ucla.edu/stat/r/dae/nbreg.htm
           newdata2 <- cbind(newdata2, predict(m1, newdata2,type =
"link", se.fit=TRUE))
           newdata2 <- within(newdata2, {
           DaysAbsent <- exp(fit)
                 LL <- exp(fit - 1.96 * se.fit)
                 UL <- exp(fit + 1.96 * se.fit)
             })
but I have an uneasy feeling about fitting standard errors from a
gaussian distribution.

In the past I've used MCMCglmm to estimate credible intervals around the
beta coefficients of a categorical variable, but now,
I have a continuous variable, conifer sapling count.

So, lots of questions!

How does one make an accurate predictive graph from a nbinom model?

Does it make sense to set up a dataframe with the means of the
covariates when the covariates'
distributions were skewed count data?

How do I use MCMCglmm correctly?

Is "se.fit" valid for the nbinom log link, or....?

my Code:
############# CONIFER SAPLINGS by pellets/ha/month (phm)  #####################
#### mod1.all <- glmmadmb(pellets ~ season * (t.con.splgs +
t.dec.trees + pctMidCov + t.BAtrees + cc) +

####                     offset(ln.days)+(1|stand/plot)+ (1|hareyr),
data=hv, family="nbinom")
predCS <- data.frame(
                     t.con.splgs = rep(seq(from = min(hv$t.con.splgs),
to = max(hv$t.con.splgs),length.out = 100), 2),

                         t.dec.trees=rep(mean(hv$t.dec.trees),200),
                         pctMidCov=rep(mean(hv$pctMidCov),200),

t.BAtrees=rep(mean(hv$t.BAtrees),200),cc=rep(mean(hv$cc),200),

                         ln.days=rep(log(30.25),200),
                         season = factor(rep(1:2, each = 100), levels
= 1:2, labels =levels(hv$season)))

predCS <- cbind(predCS, predict(mod1.all, predCS, type = "link", se.fit=TRUE))


predCS <- within(predCS, {
  pellets <- exp(fit)
  LL <- exp(fit - 1.96 * se.fit)   # these seem wrong !!
  UL <- exp(fit + 1.96 * se.fit)   # ??
})

head(predCS)           #  check to see fits

  t.con.splgs t.dec.trees pctMidCov t.BAtrees       cc  ln.days season
      fit    se.fit       UL       LL   pellets
1   0.0000000    3.493585  40.80654  4.117168 82.84937 3.409496    smr
-1.051665 0.5187013 6437.304 842.6536 0.3493557

2   0.7142493    3.493585  40.80654  4.117168 82.84937 3.409496    smr
-1.046329 0.5187194 6471.975 847.1320 0.3512248
3   1.4284985    3.493585  40.80654  4.117168 82.84937 3.409496    smr
-1.040993 0.5187633 6507.163 851.5911 0.3531040

4   2.1427478    3.493585  40.80654  4.117168 82.84937 3.409496    smr
-1.035657 0.5188331 6542.872 856.0304 0.3549932
5   2.8569971    3.493585  40.80654  4.117168 82.84937 3.409496    smr
-1.030321 0.5189286 6579.111 860.4492 0.3568926

6   3.5712464    3.493585  40.80654  4.117168 82.84937 3.409496    smr
-1.024984 0.5190500 6615.885 864.8472 0.3588020
  con.splgs.ha      phm
1      0.00000 2329.038
2      5.10152 2341.499
3     20.40608 2354.027

4     45.91368 2366.622
5     81.62432 2379.284
6    127.53801 2392.014

############# pellets/ha/month (phm) scale up
#############  and scale up saplings to per ha from 0.1ha plot level
predCS <- within(predCS, {

  pellets <- exp(fit)
  phm <- pellets/1.5*10000
   LL <- (exp(fit - 1.96 * se.fit))/1.5*10000
   UL <- (exp(fit + 1.96 * se.fit))/1.5*10000
   con.splgs.ha <- 10*(t.con.splgs^2)          # backtransform square
root and scale up

})
CS <- ggplot(predCS, aes(con.splgs.ha, phm)) +

  theme_bw() +   #eliminates background, gridlines, but not border
  theme(
    plot.background = element_blank()
   ,panel.grid.major = element_blank()

   ,panel.grid.minor = element_blank()
#   ,panel.border = element_blank()
   ,panel.background = element_blank()
  ) +

  geom_ribbon(aes(ymin = LL, ymax = UL, fill = season),alpha = .25) +
  #geom_line(aes(colour = season),size = 1) +

  geom_smooth(aes(colour=season, linetype=season),size = 1, se = F) +
  scale_colour_manual(values=c("wtr"= 4, "smr" = 3)) +
  scale_x_continuous(limits = c(0, 25000)) +
  scale_y_continuous(limits = c(0, 12000)) +

  scale_fill_brewer(palette="Accent") +

  theme(legend.position = 'none') +         #get rid of the legend
  labs(x = "Conifer Sapling Count per ha", y = "Predicted Pellets/ha/month")

print(CS)

Thanks for any ideas.
Sheryn




On Tue, Dec 10, 2013 at 11:02 PM, Ben Bolker <bbolker@gmail.com> wrote:

> On 13-12-10 08:59 PM, Sheryn Olson wrote:
> > Hello everyone,
> >
> > In this model:
> >
> > strctr.all <- glmmadmb(pellets ~ season * (totl.splgs + totl.trees +
> >                                        pctMidCov + lc +
> >                                        cc + BAsplgs + BAtrees) +
> >                                        offset(ln.days)+(1|plot/stand)+
> > (1|year),
> >                                        data=hv, family="nbinom")
> >
> > Is (1|plot/stand) correct to mean group=stand ?
>
>   There are two spatial grouping levels here -- it sounds like you
> should use stand/plot (i.e., stand and plot within stand).
> >
> > Is the repeated measurement correctly addressed with (1|year) ?
> >
> > There are 3 repeated measures of summer (2010,2011,2012); 3 of winter
> > (years 2011,2012,2013) but the 3 measures per season are a "nuisance"
> > source of random variation.   The Season variable Is the fixed variable
> of
> > interest.
>
>   You should probably use Year+(1|Season:Year), i.e. Season and Year are
> fixed (Year fixed because there aren't enough years to estimate
> among-year-variance reliably).
>
> Each plot is sampled multiple times, right?
>
>
> > Data structure: 2,184 plot records of pellet counts - the response
> variable.
> > In 26 forest stands were subsampled with 20 summer and 10 winter
> > Plots/stand for vegetation metrics. Stands are replicates.
> >
> > Thanks much,
> > Sheryn
> >
>
>
-- 
Sheryn Olson
Department of Wildlife Ecology
The University of Maine

	[[alternative HTML version deleted]]

