thr3ads.net - R help - [R] Analysing insecticide biossays using lmer [Jun 2011]

If this information is useful, please help other people find it:
Share via:

pierricklabbe

2011-Jun-30 13:56 UTC

[R] Analysing insecticide biossays using lmer

Hi all,

Here is my problem: I performed bioassays using a unique insecticide on 9
different genotypes and got their mortality depending on the dose of
insecticide used.
Now, I want to see wether some genotypes are different or not in their
responses to insecticide.
My problem is that I have up to four replicates for some genotypes, but only
one for other... Due to this unbalanced design, I thought of using lmer
procedure (lme4).

Here are the data, the model and the R output:
> dataf    geno  ml bouteille  dose   logdose tot mort date
1   D1D1 0.4      0.10 4e-04 -3.397940  32    0    1
2   D1D1 0.4      0.10 4e-04 -3.397940  83    4    2
3   D1D1 0.5      0.10 5e-04 -3.301030 120   17    2
4   D1D1 0.5      0.10 5e-04 -3.301030 114    2    4
5   D1D1 0.6      0.10 6e-04 -3.221849  34   10    1
6   D1D1 0.6      0.10 6e-04 -3.221849  72   20    2
7   D1D1 0.7      0.10 7e-04 -3.154902 120   64    2
8   D1D1 0.7      0.10 7e-04 -3.154902 116    7    4
9   D1D1 0.8      0.10 8e-04 -3.096910  35   15    1
10  D1D1 0.8      0.10 8e-04 -3.096910  78   47    2
11  D1D1 0.8      0.10 8e-04 -3.096910 120   73    3
12  D1D1 0.9      0.10 9e-04 -3.045757  32   19    1
13  D1D1 0.9      0.10 9e-04 -3.045757  87   57    2
14  D1D1 0.9      0.10 9e-04 -3.045757 119   78    3
15  D1D1 0.9      0.10 9e-04 -3.045757 116   14    4
16  D1D1 1.0      0.10 1e-03 -3.000000  33   26    1
17  D1D1 1.0      0.10 1e-03 -3.000000  83   67    2
18  D1D1 1.0      0.10 1e-03 -3.000000 120   77    3
19  D1D1 1.0      0.10 1e-03 -3.000000 120   21    4
20  D1D1 0.3      1.00 3e-03 -2.522879  38   38    1
21  D1D1 0.3      1.00 3e-03 -2.522879  94   94    2
22  D1D1 0.3      1.00 3e-03 -2.522879 120  109    3
23  D1D1 0.5      1.00 5e-03 -2.301030  40   40    1
24  D1D1 0.5      1.00 5e-03 -2.301030  94   94    2
25  D1D1 0.5      1.00 5e-03 -2.301030 120  117    3
26  D1D1 0.6      1.00 6e-03 -2.221849 120  120    2
27  D1D1 0.7      1.00 7e-03 -2.154902  40   40    1
28  D1D1 0.7      1.00 7e-03 -2.154902  94   94    2
29   D1R 0.6      0.10 6e-04 -3.221849  87   25    1
30   D1R 0.8      0.10 8e-04 -3.096910  88   40    1
31   D1R 1.0      0.10 1e-03 -3.000000  88   39    1
32   D1R 0.3      1.00 3e-03 -2.522879  93   59    1
33   D1R 0.4      1.00 4e-03 -2.397940  96   85    1
34   D1R 0.5      1.00 5e-03 -2.301030  96   91    1
35   D1R 0.7      1.00 7e-03 -2.154902  95   91    1
36   D1R 0.9      1.00 9e-03 -2.045757  97   97    1
37   D1S 0.3      0.10 3e-04 -3.522879 117    4    1
38   D1S 0.5      0.10 5e-04 -3.301030 118   37    1
39   D1S 0.6      0.10 6e-04 -3.221849 119   76    1
40   D1S 0.7      0.10 7e-04 -3.154902 118   98    1
41   D1S 0.9      0.10 9e-04 -3.045757 120  116    1
42   D1S 1.0      0.10 1e-03 -3.000000 120  120    1
43  D3D3 0.5      0.10 5e-04 -3.301030  58    6    2
44  D3D3 0.5      0.10 5e-04 -3.301030  40    4    3
45  D3D3 0.7      0.10 7e-04 -3.154902  57   12    2
46  D3D3 0.7      0.10 7e-04 -3.154902  40    9    3
47  D3D3 0.8      0.10 8e-04 -3.096910  87   13    1
48  D3D3 0.8      0.10 8e-04 -3.096910  59   12    2
49  D3D3 0.9      0.10 9e-04 -3.045757  97   35    1
50  D3D3 0.9      0.10 9e-04 -3.045757  56   13    2
51  D3D3 0.9      0.10 9e-04 -3.045757  39   11    3
52  D3D3 1.0      0.10 1e-03 -3.000000  96   49    1
53  D3D3 1.0      0.10 1e-03 -3.000000  60   20    2
54  D3D3 1.0      0.10 1e-03 -3.000000  35    6    3
55  D3D3 0.2      1.00 2e-03 -2.698970  40   31    1
56  D3D3 0.3      1.00 3e-03 -2.522879 100   88    1
57  D3D3 0.3      1.00 3e-03 -2.522879  60   50    2
58  D3D3 0.3      1.00 3e-03 -2.522879  40   37    3
59  D3D3 0.5      1.00 5e-03 -2.301030  99   93    1
60  D3D3 0.5      1.00 5e-03 -2.301030  60   59    2
61  D3D3 0.5      1.00 5e-03 -2.301030  40   39    3
62  D3D3 0.6      1.00 6e-03 -2.221849 100  100    1
63  D3D3 0.6      1.00 6e-03 -2.221849  60   60    2
64  D3D3 0.6      1.00 6e-03 -2.221849  40   40    3
65   D3R 0.6      0.10 6e-04 -3.221849  87    7    1
66   D3R 0.8      0.10 8e-04 -3.096910  84   30    1
67   D3R 1.0      0.10 1e-03 -3.000000  87   36    1
68   D3R 0.3      1.00 3e-03 -2.522879  92   73    1
69   D3R 0.4      1.00 4e-03 -2.397940 100   91    1
70   D3R 0.5      1.00 5e-03 -2.301030  99   95    1
71   D3R 0.7      1.00 7e-03 -2.154902  94   93    1
72   D3R 0.8      1.00 8e-03 -2.096910  98   98    1
73   D3S 1.0      0.01 1e-04 -4.000000  30    0    1
74   D3S 1.0      0.01 1e-04 -4.000000  60    0    2
75   D3S 1.0      0.01 1e-04 -4.000000  61    0    3
76   D3S 0.3      0.10 3e-04 -3.522879  28    1    1
77   D3S 0.3      0.10 3e-04 -3.522879  66    1    3
78   D3S 0.4      0.10 4e-04 -3.397940  66    2    3
79   D3S 0.5      0.10 5e-04 -3.301030  40    9    1
80   D3S 0.5      0.10 5e-04 -3.301030  60   57    2
81   D3S 0.5      0.10 5e-04 -3.301030  67    8    3
82   D3S 0.6      0.10 6e-04 -3.221849  42   12    1
83   D3S 0.6      0.10 6e-04 -3.221849  60   59    2
84   D3S 0.6      0.10 6e-04 -3.221849  66   22    3
85   D3S 0.7      0.10 7e-04 -3.154902  38   25    1
86   D3S 0.7      0.10 7e-04 -3.154902  60   60    2
87   D3S 0.7      0.10 7e-04 -3.154902  70   36    3
88   D3S 0.8      0.10 8e-04 -3.096910  72   51    3
89   D3S 0.9      0.10 9e-04 -3.045757  45   43    1
90   D3S 1.0      0.10 1e-03 -3.000000  48   48    1
91    RR 0.2      1.00 2e-03 -2.698970  39    1    1
92    RR 0.2      1.00 2e-03 -2.698970  60    0    2
93    RR 0.4      1.00 4e-03 -2.397940  39    3    1
94    RR 0.4      1.00 4e-03 -2.397940  60    5    2
95    RR 0.5      1.00 5e-03 -2.301030  60    9    2
96    RR 0.6      1.00 6e-03 -2.221849  39   13    1
97    RR 0.6      1.00 6e-03 -2.221849  60   22    2
98    RR 0.7      1.00 7e-03 -2.154902  60   42    2
99    RR 0.8      1.00 8e-03 -2.096910  40   38    1
100   RR 0.8      1.00 8e-03 -2.096910  60   39    2
101   RR 0.2     10.00 2e-02 -1.698970  40   40    1
102   RR 0.2     10.00 2e-02 -1.698970  60   60    2
103   RS 0.5      0.10 5e-04 -3.301030  38    0    1
104   RS 0.6      0.10 6e-04 -3.221849  40    2    1
105   RS 0.8      0.10 8e-04 -3.096910  40    1    1
106   RS 0.9      0.10 9e-04 -3.045757 114   21    1
107   RS 1.0      0.10 1e-03 -3.000000  39    1    1
108   RS 0.2      1.00 2e-03 -2.698970 120   82    1
109   RS 0.3      1.00 3e-03 -2.522879  38   33    1
110   RS 0.4      1.00 4e-03 -2.397940 120  110    1
111   RS 0.6      1.00 6e-03 -2.221849 120  119    1
112   RS 0.8      1.00 8e-03 -2.096910  40   40    1
113   SS 0.8      0.01 8e-05 -4.096910  72    0    2
114   SS 0.9      0.01 9e-05 -4.045757  37    9    1
115   SS 0.9      0.01 9e-05 -4.045757  75    1    2
116   SS 1.0      0.01 1e-04 -4.000000  39    5    1
117   SS 1.0      0.01 1e-04 -4.000000  73    1    2
118   SS 0.2      0.10 2e-04 -3.698970  40   21    1
119   SS 0.2      0.10 2e-04 -3.698970  77   12    2
120   SS 0.3      0.10 3e-04 -3.522879  38   28    1
121   SS 0.4      0.10 4e-04 -3.397940  40   30    1
122   SS 0.5      0.10 5e-04 -3.301030  39   37    1
123   SS 0.5      0.10 5e-04 -3.301030  41   36    2
124   SS 0.6      0.10 6e-04 -3.221849  40   40    1
125   SS 0.6      0.10 6e-04 -3.221849  40   39    2> attach(dataf)
> y<-cbind(mort,tot-mort)
>
model1<-lmer(y~logdose*geno+(1|geno/logdose)-1,family="quasibinomial",data=dataf)Erreur : length(f1) == length(f2) is not TRUE
De plus : Warning messages:
1: In logdose:geno :
  l'expression num?rique a 125 ?l?ments : seul le premier est utilis?
2: In logdose:geno :
  l'expression num?rique a 125 ?l?ments : seul le premier est utilis?

Enter a frame number, or 0 to exit   

 1: lmer(y ~ logdose * geno + (1 | geno/logdose) - 1, family
"quasibinomial", data = dataf)
 2: eval.parent(mc)
 3: eval(expr, p)
 4: eval(expr, envir, enclos)
 5: glmer(formula = y ~ logdose * geno + (1 | geno/logdose) - 1, data dataf,
family = "quasibinomi
 6: lmerFactorList(formula, fr, 0, 0)
 7: sapply(seq_along(fl)[-1], function(i) isNested(fl[[i - 1]], fl[[i]]))
 8: lapply(X, FUN, ...)
 9: FUN(2[[1]], ...)
10: isNested(fl[[i - 1]], fl[[i]])
11: stopifnot(length(f1) == length(f2))

S?lection : summary (model1)
Enter an item from the menu, or 0 to exit
S?lection : 10
Called from: eval(expr, envir, enclos)
Browse[1]> f1
  [1] D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1
D1D1 D1D1 D1D1 D1D1 D1D1
 [20] D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1D1 D1R  D1R  D1R  D1R  D1R 
D1R  D1R  D1R  D1S  D1S 
 [39] D1S  D1S  D1S  D1S  D3D3 D3D3 D3D3 D3D3 D3D3 D3D3 D3D3 D3D3 D3D3 D3D3
D3D3 D3D3 D3D3 D3D3 D3D3
 [58] D3D3 D3D3 D3D3 D3D3 D3D3 D3D3 D3D3 D3R  D3R  D3R  D3R  D3R  D3R  D3R 
D3R  D3S  D3S  D3S  D3S 
 [77] D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S  D3S 
RR   RR   RR   RR   RR  
 [96] RR   RR   RR   RR   RR   RR   RR   RS   RS   RS   RS   RS   RS   RS  
RS   RS   RS   SS   SS  
[115] SS   SS   SS   SS   SS   SS   SS   SS   SS   SS   SS  
Levels: D1D1 D1R D1S D3D3 D3R D3S RR RS SS
Browse[1]> f2
[1] -3.397940009 -2.397940009 -1.397940009 -0.397940009 0.602059991 
Levels: -3.397940009 -2.397940009 -1.397940009 -0.397940009 0.602059991

As you can see, the model does not work because it does not consider all the
logdose values...
Could someone help me in finding what is wrong or in suggesting maybe
another path to follow for my statistical analysis?

Cheers

Pierrick


--
View this message in context:
http://r.789695.n4.nabble.com/Analysing-insecticide-biossays-using-lmer-tp3635724p3635724.html
Sent from the R help mailing list archive at Nabble.com.

Ben Bolker

2011-Jun-30 18:48 UTC

head link

[R] Analysing insecticide biossays using lmer

pierricklabbe <pierrick.labbe <at> univ-montp2.fr> writes:
> 
> Hi all,
> 
> Here is my problem: I performed bioassays using a unique insecticide on 9
> different genotypes and got their mortality depending on the dose of
> insecticide used.
> Now, I want to see wether some genotypes are different or not in their
> responses to insecticide.
> My problem is that I have up to four replicates for some genotypes, but
only
> one for other... Due to this unbalanced design, I thought of using lmer
> procedure (lme4).
> 
  There are a few issues here.

(1) please send follow-ups to r-sig-mixed-models at r-project.org
(2) quasi- families are not available in recent versions of lme4:
you must have a pretty old version.  The current advice is to use
an observation-level random effect to capture overdispersion in GLMMs.
(3) the variables used in defining groups (i.e. the ones found to the
right of the | in the random-effects part of the formula) must be
factors.  If you want a random effect of the continuous dose variable
you need to put it on the *left* of the bar (this is the proximate cause
of your error message).  You could alternatively make logdose into
a factor, but with 18 unique values I don't think that's a good idea.

dataf <- read.table("insect.dat",header=TRUE)
library(lme4)

dataf$obs <- seq(nrow(dataf))

model1 <- lmer(cbind(mort,tot-mort)~logdose*geno-1+(logdose|geno)+(1|obs),
               family=binomial,data=dataf)

model0 <- update(model1,.~.-(1|obs))

AIC(model0,model1)

By the way, you should always look at your data ...

library(ggplot2)
library(mgcv)
ggplot(dataf,aes(x=logdose,y=mort/tot,colour=geno))+
  geom_point()+geom_smooth(method="gam",family=binomial,
     aes(weight=tot))

Maybe Matching Threads

Search for more apparently analagous threads

R help - Jun 2011 - Analysing insecticide biossays using lmer

[R] Analysing insecticide biossays using lmer

[R] Analysing insecticide biossays using lmer

Maybe Matching Threads