[R-sig-ME] glmer question family=Gamma(link="identity") error. Re: lmer formula syntax?

[Colin Wahl]

Generalized mixed model geeks,

I just recently asked a question to this list about lmer formula syntax and
was fortunate to be flattered by your helpful responses (Re: lmer formula
syntax). I now have a new question regarding the specification of my error
distribution.

My previous study description: I am performing an ecological study of stream
health using %ept (An aggregate relative abundance variable for
invertebrates) as my response variable. I am interested in how invertebrate
counts vary between different watershed types (forested, cultivated,
developed, grassland) and reach types (forested, non forested). There are a
total of 12 streams, 12 watersheds (wsh) and 24 riparia(rip). I am using OSx


> summary(ept[,2:5])
 wsh        stream       rip         ept
 c:18        BA     : 6   F:36      Min.   :0.000000
 d:18        D1     : 6   N:36     1st Qu.:0.005797
 f:24         D2     : 6                Median :0.029630
 g:12        DC     : 6               Mean   :0.117869
                FE     : 6               3rd Qu.:0.221027
                JO     : 6               Max.   :0.527090
               (Other):36

The complete data set is detailed below.
stream is the only random variable and it is nested (not implicitly) in
wshed. My statistical model is as follows:
ept(ijkl) = μ + wshedi  + stream(i)j + ripariak + wshed*ripariaik +
stream*riparia(i)jk + ε(ijk)l

My response variable (%ept) is the proportion of invertebrates sampled that
belong to sensitive orders (ephemeroptera, plecoptera, tricoptera). This is
a binomial response ("ept" or "not ept"), but I am using a proportion
instead of discrete counts. The error and variance structures should be
better modeled with a binomial family (or Gamma for continuous [or
proportional] data?) I tried this using the following formula with a Gamma
distribution and an identity link:

>model<-glmer(ept ~ wsh*rip + (1|stream) + (1|stream:rip), data=ept,
family=Gamma(link="identity"))

...and received this error:
Error in eval(expr, envir, enclos) :
  non-positive values not allowed for the gamma family

Then I cleverly tried to trick it by adding 0.1 to "ept" and was told in a
strange language that my clever tricks are just feeble ditherings:
Error in mer_finalize(ans) :
  mu[i] must be positive: mu = -0.0120754, i = 9

Does this mean that the proposed model is estimating negative values? None
of my response values are negative. There must be some way to correctly
model my data with an appropriate binomial-like distribution. The error
structure would be vastly improved.

Do any of you helpful (g)lmeronauts have any suggestions on how I can do
this?

It was suggested to me that there may be some way to specify a binomial
distribution in such a way that proportions are an acceptable form of
response (because this is a necessary part of the model estimation anyway,
right?)

Thank you!
Colin

Detailed data

> ept[,2:5]
   wsh stream rip     ept
1    c     DC   F 0.09770
2    c     DC   F 0.02419
3    c     DC   F 0.03704
4    c     JO   F 0.02128
5    c     JO   F 0.02478
6    c     JO   F 0.02358
7    c     SD   F 0.00000
8    c     SD   F 0.00221
9    c     SD   F 0.00000
10   c     DC   N 0.01794
11   c     DC   N 0.01235
12   c     DC   N 0.01717
13   c     JO   N 0.05556
14   c     JO   N 0.02000
15   c     JO   N 0.00000
16   c     SD   N 0.00591
17   c     SD   N 0.02941
18   c     SD   N 0.02985
19   d     FE   F 0.00293
20   d     FE   F 0.00218
21   d     FE   F 0.00749
22   d     PA   F 0.00000
23   d     PA   F 0.00000
24   d     PA   F 0.00000
25   d     SC   F 0.00073
26   d     SC   F 0.00673
27   d     SC   F 0.00173
28   d     FE   N 0.00000
29   d     FE   N 0.00000
30   d     FE   N 0.00000
31   d     PA   N 0.00781
32   d     PA   N 0.00182
33   d     PA   N 0.00684
34   d     SC   N 0.00000
35   d     SC   N 0.03604
36   d     SC   N 0.00000
37   f     BA   F 0.22841
38   f     BA   F 0.30508
39   f     BA   F 0.16864
40   f     D1   F 0.45753
41   f     D1   F 0.41422
42   f     D1   F 0.46531
43   f     D2   F 0.28852
44   f     D2   F 0.22123
45   f     D2   F 0.25141
46   f     SP   F 0.28907
47   f     SP   F 0.14908
48   f     SP   F 0.26863
49   f     BA   N 0.13995
50   f     BA   N 0.34842
51   f     BA   N 0.39713
52   f     D1   N 0.26800
53   f     D1   N 0.23448
54   f     D1   N 0.19345
55   f     D2   N 0.32421
56   f     D2   N 0.22096
57   f     D2   N 0.21422
58   f     SP   N 0.03133
59   f     SP   N 0.02703
60   f     SP   N 0.04032
61   g     MC   F 0.19917
62   g     MC   F 0.16966
63   g     MC   F 0.15318
64   g     SI   F 0.33650
65   g     SI   F 0.52709
66   g     SI   F 0.47359
67   g     MC   N 0.01029
68   g     MC   N 0.00546
69   g     MC   N 0.00816
70   g     SI   N 0.02693
71   g     SI   N 0.09756
72   g     SI   N 0.03907

> summary(ept[,2:5])
 wsh        stream       rip         ept
 c:18        BA     : 6   F:36      Min.   :0.000000
 d:18        D1     : 6   N:36     1st Qu.:0.005797
 f:24         D2     : 6                Median :0.029630
 g:12        DC     : 6               Mean   :0.117869
                FE     : 6               3rd Qu.:0.221027
                JO     : 6               Max.   :0.527090
               (Other):36

> str(ept[,2:5])
'data.frame':    72 obs. of  4 variables:
 $ wsh   : Factor w/ 4 levels "c","d","f","g": 1 1 1 1 1 1 1 1 1 1 ...
 $ stream: Factor w/ 12 levels "BA","D1","D2",..: 4 4 4 6 6 6 10 10 10 4 ...
 $ rip   : Factor w/ 2 levels "F","N": 1 1 1 1 1 1 1 1 1 2 ...
 $ ept   : num  0.0977 0.0242 0.037 0.0213 0.0248 ...

-- 
Colin Wahl
Department of Biology
Western Washington University
Bellingham WA, 98225
ph: 360-391-9881