Plotting drm in DRC package - multiple model plot

Question

1. Curves in DRC package not fitting correctly on plot
2. Broken = TRUE graphical parameter appearing outside the plot
3. Applying type = "confidence" yields an error

Hi all,

I am having trouble plotting a drm in the R package DRC.

My data comprises of three variables; Concentration (numeric vector of concentrations of compound used in assays in mg/L), Herbicide (Character vector of type of compound), Inhib (numeric vector of inhibition per test treatment).

The dose response curve works fine, and individually the plots and models work.

The problem is that the curve on two of the compounds continues above the maximum inhibition for that herbicide.

I rescaled the y axis and it does not fix the problem:

I tried re-organising the dataframe so that a different herbicide was at the top instead (i.e D or E instead of A or B) and the continuing ill-fitting curve seems to effect the first two herbicides whatever they are.

Thanks in advance!

My code is seen below and the error (Q 3) is the last line:

multi.m1 <- drm(mydata$Inhib~mydata$Concentration, data = mydata, Herbicide, fct = LL.4())
plot(multi.m1) ##try limit y axis, change x axis for broken

plot(multi.m1, col = TRUE, pch = 20:25, xlab = "Concentration", ylab = "Inhibition", broken = TRUE)

plot(multi.m1, ylim = c(-5, 102), broken = TRUE, add = TRUE, 
 type = "confidence") #why does x axis break appear wonky? -> ask the oracle

Error in parmMat[, groupLevels, drop = FALSE] : subscript out of bounds

EDIT, mydata:

 Concentration Herbicide       Inhib
1      2.2375e-02         C  84.1171273
2      2.2375e-02         C  83.2708908
3      2.2375e-02         C  80.7519653
4      8.9500e-03         C  73.7143264
5      8.9500e-03         C  76.9312856
6      8.9500e-03         C  69.5508871
7      3.5800e-03         C  58.9470598
8      3.5800e-03         C  60.3750236
9      3.5800e-03         C  58.1370479
10     1.4320e-03         C  30.2296338
11     1.4320e-03         C  26.6788108
12     1.4320e-03         C  45.5096997
13     5.7300e-04         C   9.0653245
14     5.7300e-04         C  -3.1497620
15     5.7300e-04         C   2.0924363
16     2.2375e-02         C  62.8053281
17     2.2375e-02         C  62.8053281
18     2.2375e-02         C  69.8827947
19     8.9500e-03         C  59.6812753
20     8.9500e-03         C  59.6812753
21     8.9500e-03         C  59.6812753
22     3.5800e-03         C  49.1470273
23     3.5800e-03         C  50.3237083
24     3.5800e-03         C  45.8087647
25     1.4320e-03         C  27.2677293
26     1.4320e-03         C  22.1870253
27     1.4320e-03         C  30.0683330
28     5.7300e-04         C  15.5899507
29     5.7300e-04         C   6.5319619
30     5.7300e-04         C   5.7281208
31     5.6200e+02         D  63.9535852
32     5.6200e+02         D  61.6934485
33     5.6200e+02         D  63.9535852
34     2.2500e+02         D  59.5337963
35     2.2500e+02         D  57.4660724
36     2.2500e+02         D  66.3240151
37     9.0000e+01         D  41.9942602
38     9.0000e+01         D  45.0317602
39     9.0000e+01         D  46.6248268
40     3.6000e+01         D  21.8067846
41     3.6000e+01         D  21.8067846
42     3.6000e+01         D  24.8129501
43     1.4400e+01         D  13.7654019
44     1.4400e+01         D  10.5784119
45     1.4400e+01         D  12.9488419
46     5.6200e+02         D 101.8769351
47     5.6200e+02         D 102.1410258
48     5.6200e+02         D  93.7975216
49     2.2500e+02         D  93.5444549
50     2.2500e+02         D  92.1264940
51     2.2500e+02         D  92.4857338
52     9.0000e+01         D  75.5170610
53     9.0000e+01         D  71.5852452
54     9.0000e+01         D  73.0984175
55     3.6000e+01         D  53.5097851
56     3.6000e+01         D  54.2592274
57     3.6000e+01         D  60.2304326
58     1.4400e+01         D  34.6248616
59     1.4400e+01         D  28.6424632
60     1.4400e+01         D  30.8926163
61     4.0000e-02         E  71.3055968
62     4.0000e-02         E  67.3600153
63     4.0000e-02         E  63.7579815
64     2.8571e-02         E  57.3765637
65     2.8571e-02         E  49.3390342
66     2.8571e-02         E  60.4444317
67     2.0408e-02         E  44.7345962
68     2.0408e-02         E  43.6592819
69     2.0408e-02         E  36.8254304
70     1.4577e-02         E  32.5538843
71     1.4577e-02         E  25.1417758
72     1.4577e-02         E  32.5538843
73     1.0472e-02         E  13.4025320
74     1.0472e-02         E  23.8062641
75     1.0472e-02         E  17.7074226
76     4.0000e-02         E  91.6851862
77     4.0000e-02         E  75.5919138
78     4.0000e-02         E  91.5043781
79     2.8571e-02         E  66.4345965
80     2.8571e-02         E  73.3145119
81     2.8571e-02         E  62.6140016
82     2.0408e-02         E  52.0485485
83     2.0408e-02         E  56.1078884
84     2.0408e-02         E  49.0660567
85     1.4577e-02         E  50.9756413
86     1.4577e-02         E  47.8435190
87     1.4577e-02         E  39.5595697
88     1.0472e-02         E  36.4083394
89     1.0472e-02         E  27.2520151
90     1.0472e-02         E  36.8735853
91     2.5000e-01         A  59.6812753
92     2.5000e-01         A  66.1900088
93     2.5000e-01         A  59.6812753
94     1.0000e-01         A  41.7461116
95     1.0000e-01         A  48.0034967
96     1.0000e-01         A  50.3237083
97     4.0000e-02         A  33.8648039
98     4.0000e-02         A  16.1005569
99     4.0000e-02         A  23.3988362
100    1.6000e-02         A  13.6094285
101    1.6000e-02         A   5.7281208
102    1.6000e-02         A   9.4777989
103    6.4000e-03         A   2.6630689
104    6.4000e-03         A   0.1587632
105    6.4000e-03         A  -0.1867083
106    2.5000e-01         A  75.2770560
107    2.5000e-01         A  74.4016723
108    2.5000e-01         A  74.0126249
109    1.0000e-01         A  51.1490611
110    1.0000e-01         A  54.9412310
111    1.0000e-01         A  52.9289775
112    4.0000e-02         A  36.6921613
113    4.0000e-02         A  20.3596646
114    4.0000e-02         A  28.7826960
115    1.6000e-02         A   7.9617212
116    1.6000e-02         A   3.5240504
117    1.6000e-02         A   1.4995895
118    6.4000e-03         A  -4.7738737
119    6.4000e-03         A -11.5182880
120    6.4000e-03         A  -6.8625243
121    2.1020e-02         B  65.1307269
122    2.1020e-02         B  63.7534772
123    2.1020e-02         B  69.6048270
124    1.4500e-02         B  39.1407346
125    1.4500e-02         B  37.1314501
126    1.4500e-02         B  36.4855432
127    1.0000e-02         B  31.1355743
128    1.0000e-02         B  25.5101867
129    1.0000e-02         B  35.2270466
130    6.8965e-03         B  20.2249022
131    6.8965e-03         B  11.1906726
132    6.8965e-03         B  21.9030270
133    4.7500e-03         B   9.9321760
134    4.7500e-03         B  10.8719563
135    4.7500e-03         B  11.1906726
136    2.1020e-02         B  83.0444344
137    2.1020e-02         B  80.0717309
138    2.1020e-02         B  83.0444344
139    1.4500e-02         B  64.4239326
140    1.4500e-02         B  68.1575575
141    1.4500e-02         B  69.0731273
142    1.0000e-02         B  42.8229958
143    1.0000e-02         B  48.0750454
144    1.0000e-02         B  58.4858694
145    6.8965e-03         B  48.3924985
146    6.8965e-03         B  47.2534892
147    6.8965e-03         B  32.4004160
148    4.7500e-03         B  25.4700950
149    4.7500e-03         B  31.1661066
150    4.7500e-03         B  28.5676908

Answer 1

Fixing the upper-limit at 100 within the model may exacerbate issues with over-fitting and may also be inconsistent with the biology that is being studied. There are two other options regarding the upper-limit which may be more appropriate.

The first is to set a maximum limit for the upper-bound using upperl;

multi.m2 <- drm(Inhib~Concentration, data = mydata, Herbicide, 
            fct = LL.4(),
            upperl=c(NA,NA,100,NA))

The second option is to use pmodels to direct drm to calculate the upper limit from the pooled data for all herbicides.

multi.m3 <- drm(Inhib~Concentration, data = mydata, Herbicide, 
            fct = LL.4(),
            pmodels=data.fram(Herbicide,Herbicide,1,Herbicide))

Answer 2

Apparently there is no specific default of 100 for an upper limit as you expect. But it's not that difficult to impose one after you follow through the trail of documentation in:

?drm
?LL.4
?ryegrass  # where you can see the naming of parameters of LL.4 parameters

So setting only the third parameter of the limits on LL.4 to your expected 100, we get:

multi.m1 <- drm(Inhib~Concentration, data = mydata[c(3,1,2)], Herbicide, 
   fct = LL.4(fixed=c(NA, NA, 100, NA), 
              names = c("Slope", "Lower Limit", "Upper Limit", "ED50")))

plot(multi.m1, col = TRUE, pch = 20:25, xlab = "Concentration", ylab = "Inhibition")

If you need to also impose a zero on the lower limit the way forward is clear. The Oracle has spoken.