Philippe Grosjean
2001-May-01 11:05 UTC
[R] SSfpl self-start sometimes fails... workaround proposed
Hello,
nls library provides 6 self-starting models, among them: SSfp, a four
parameters logistic function. Its self-starting procedure involves several
steps. One of these steps is:
pars <- as.vector(coef(nls(y ~ cbind(1, 1/(1 + exp((xmid - x)/exp(lscal)))),
data = xydata, start = list(lscal = 0), algorithm = "plinear")))
which assumes an initial value of lscal equal to 0. If lscal is very
different to 0, the evaluation could fail (singular gradient,...), as it is
the case with the dataset provided hereunder (see end of this message).
As a workaround, I propose a modified self-start function for SSfpl that
remedy to this:
---------------------------------------
fpl <- function (input, A, B, xmid, scal)
{
.expr1 <- B - A
.expr2 <- xmid - input
.expr4 <- exp((.expr2/scal))
.expr5 <- 1 + .expr4
.expr8 <- 1/.expr5
.expr13 <- .expr5^2
.value <- A + (.expr1/.expr5)
.actualArgs <- as.list(match.call()[c("A", "B",
"xmid", "scal")])
if (all(unlist(lapply(.actualArgs, is.name)))) {
.grad <- array(0, c(length(.value), 4), list(NULL, c("A",
"B",
"xmid", "scal")))
.grad[, "A"] <- 1 - .expr8
.grad[, "B"] <- .expr8
.grad[, "xmid"] <- -((.expr1 * (.expr4 *
(1/scal)))/.expr13)
.grad[, "scal"] <- (.expr1 * (.expr4 *
(.expr2/(scal^2))))/.expr13
dimnames(.grad) <- list(NULL, .actualArgs)
attr(.value, "gradient") <- .grad
}
.value
}
# Self-starting function for fpl
fpl.ival <- function (mCall, data, LHS)
{
xy <- sortedXyData(mCall[["input"]], LHS, data)
if (nrow(xy) < 5) {
stop("Too few distinct input values to fit a four-parameter
logistic")
}
xydata <- c(as.list(xy), xmid = NLSstClosestX(xy,
mean(range(xy[["y"]]))))
xydata <- as.list(xydata)
options(show.error.messages = FALSE)
# Sometimes, the following evaluation gives an error (singular gradient,
...). Try another way to estimate initial parameters
pars <- try(as.vector(coef(nls(y ~ cbind(1, 1/(1 + exp((xmid -
x)/exp(lscal)))), data = xydata, start = list(lscal = 0), algorithm
"plinear"))))
if (!is.null(class(pars)) && class(pars) == "try-error") {
cat("Using second alternative for initial parameters
estimation\n")
options(show.error.messages = TRUE)
pars <- as.vector(coef(nls(y ~ SSlogis(x, Asym, xmid, scal),
data=xydata)))
xydata$xmid <- pars[2]
pars[1] <- log(pars[3])
} else { options(show.error.messages = TRUE)}
pars <- as.vector(coef(nls(y ~ cbind(1, 1/(1 + exp((xmid -
x)/exp(lscal)))), data = data.frame(xy), start = list(xmid = xydata$xmid,
lscal = pars[1]), algorithm = "plinear")))
value <- c(pars[3], pars[3] + pars[4], pars[1], exp(pars[2]))
names(value) <- mCall[c("A", "B", "xmid",
"scal")]
value
}
# Self-Starting four-parameter logistic, combining fpl and fpl.ival
SSfpl <- selfStart(fpl, fpl.ival)
remove(fpl, fpl.ival)
attr(SSfpl,"pnames") <- c("A", "B",
"xmid", "scal")
---------------------------------------
Now the self-starting function succeed with the faulty dataset (if it is
loaded under "Sample", we got):> Sample.nls <- nls(y ~ SSfpl(x, A, B, xmid, scal), data=Sample, trace=T)
Using second alternative for initial parameters estimation
3912.605 : -6.838615 97.175628 1079.915732 345.051411
This method succeed also when A is very different to 0, as
in:> Sample$y <- Sample$y - 30
>Sample.nls <- nls(y ~ SSfpl(x, A, B, xmid, scal), data=Sample, trace=T)
Using second alternative for initial parameters estimation
3912.605 : -36.83861 67.17563 1079.91578 345.05135> Sample$y <- Sample$y + 60
> Sample.nls <- nls(y ~ SSfpl(x, A, B, xmid, scal), data=Sample, trace=T)
Using second alternative for initial parameters estimation
3912.605 : 23.16139 127.17563 1079.91574 345.05139
Here is the Sample dataset that causes failure in original self-start
function of SSfpl:
--------------------------------------------
x y
1 100 3.0993441
2 200 7.1703742
3 300 2.6698554
4 400 -1.1030135
5 500 3.1997802
6 600 13.4839286
7 700 17.5247783
8 800 25.3965763
9 900 32.7809024
10 1000 42.8638906
11 1100 49.7806380
12 1200 53.3720730
13 1300 57.1732873
14 1400 69.1322680
15 1500 74.2821286
16 1600 84.7223610
17 1700 87.1016435
18 1800 83.7785903
19 1900 82.6630488
20 2000 85.6699881
21 2100 93.6834595
22 2200 87.6399878
23 2300 98.0706325
24 2400 87.7090574
25 2500 87.8427111
26 100 6.1086982
27 200 -0.2703793
28 300 1.1789864
29 400 8.4005067
30 500 7.5632783
31 600 16.4558476
32 700 15.3683658
33 800 25.4356474
34 900 36.7496044
35 1000 46.2637329
36 1100 49.4245916
37 1200 49.2426203
38 1300 67.6765838
39 1400 72.1293952
40 1500 72.6686691
41 1600 73.1239807
42 1700 83.2326747
43 1800 83.4350708
44 1900 92.0119974
45 2000 88.1731913
46 2100 92.6361306
47 2200 87.6798161
48 2300 89.8388569
49 2400 99.9118565
50 2500 94.0181265
51 100 -8.5761223
52 200 -2.8684148
53 300 -5.0846235
54 400 8.0067485
55 500 15.9500212
56 600 7.8298746
57 700 14.8714996
58 800 31.1773970
59 900 25.5603584
60 1000 36.6270604
61 1100 49.1396008
62 1200 55.1484732
63 1300 55.6742251
64 1400 63.0085688
65 1500 73.7323245
66 1600 79.7760679
67 1700 78.0336433
68 1800 88.3090297
69 1900 91.4579878
70 2000 89.4620021
71 2100 86.0521263
72 2200 92.4597165
73 2300 102.3826961
74 2400 89.6810273
75 2500 104.4345010
76 100 4.1908793
77 200 1.1663101
78 300 1.2138565
79 400 4.6746695
80 500 10.6366307
81 600 13.5696400
82 700 19.4830518
83 800 24.9596299
84 900 39.2800245
85 1000 33.1233443
86 1100 45.4346350
87 1200 58.0112784
88 1300 57.4907772
89 1400 66.8782190
90 1500 76.1078100
91 1600 77.9702662
92 1700 86.4002619
93 1800 90.8551005
94 1900 93.0494288
95 2000 91.3292011
96 2100 92.8803548
97 2200 95.3348279
98 2300 93.3454901
99 2400 89.4888304
100 2500 101.6825561
101 100 -2.6463179
102 200 2.4626755
103 300 8.8313503
104 400 8.1322231
105 500 11.2579274
106 600 7.6675020
107 700 22.8678546
108 800 18.4455407
109 900 30.8940806
110 1000 41.5778433
111 1100 44.4458981
112 1200 58.2412531
113 1300 60.1528031
114 1400 66.8289217
115 1500 74.2174892
116 1600 82.1478187
117 1700 87.0719569
118 1800 86.3247909
119 1900 77.6874618
120 2000 89.9793056
121 2100 91.5725152
122 2200 86.8773791
123 2300 90.9528753
124 2400 94.8046681
125 2500 96.2915658
126 100 3.3783906
127 200 -3.7782437
128 300 -5.3800214
129 400 8.5946633
130 500 10.1502830
131 600 14.4726541
132 700 25.4276710
133 800 26.1740854
134 900 32.4392613
135 1000 36.8220085
136 1100 45.3711596
137 1200 55.1290544
138 1300 60.0386737
139 1400 72.9271852
140 1500 74.7238370
141 1600 78.0457813
142 1700 73.7611185
143 1800 90.1935669
144 1900 89.8282116
145 2000 93.3865092
146 2100 96.5737605
147 2200 91.2951522
148 2300 93.1026883
149 2400 97.9327117
150 2500 96.5161475
151 100 -3.1399066
152 200 4.2818480
153 300 2.1415215
154 400 6.9593208
155 500 5.7589779
156 600 16.9917719
157 700 20.5212153
158 800 24.9654495
159 900 41.4699654
160 1000 44.7710150
161 1100 42.1270085
162 1200 57.5482135
163 1300 57.2669363
164 1400 68.3427189
165 1500 68.6810161
166 1600 71.9176212
167 1700 81.1582222
168 1800 92.1411063
169 1900 88.6745264
170 2000 95.1665112
171 2100 91.1493740
172 2200 93.8791498
173 2300 98.9862266
174 2400 94.7733683
175 2500 101.4603052
176 100 4.0875772
177 200 -2.5260400
178 300 -2.4073486
179 400 5.3199050
180 500 10.4923312
181 600 17.7314494
182 700 14.8992854
183 800 28.1606941
184 900 31.1065127
185 1000 38.3060603
186 1100 43.7539674
187 1200 55.7475726
188 1300 52.3597949
189 1400 67.9419737
190 1500 69.4292643
191 1600 78.9012221
192 1700 83.6830419
193 1800 93.1557145
194 1900 92.4328677
195 2000 81.7399979
196 2100 95.5528195
197 2200 94.2200239
198 2300 91.2460136
199 2400 96.5153395
200 2500 99.5341194
201 100 1.0072133
202 200 1.9475437
203 300 0.8895067
204 400 7.2445463
205 500 5.8782348
206 600 10.8238589
207 700 16.0272159
208 800 19.8050738
209 900 36.4075759
210 1000 44.4547541
211 1100 43.7834169
212 1200 55.6825097
213 1300 61.4506170
214 1400 73.7791672
215 1500 68.0171640
216 1600 70.7886840
217 1700 85.1198613
218 1800 85.9437649
219 1900 91.6470669
220 2000 90.6371665
221 2100 89.6053979
222 2200 90.6202220
223 2300 95.5400639
224 2400 93.7101926
225 2500 98.8553117
226 100 -3.2451954
227 200 4.0275414
228 300 5.3651481
229 400 1.9235670
230 500 12.6634105
231 600 19.4137491
232 700 13.5260270
233 800 25.0315762
234 900 30.3619952
235 1000 37.0424986
236 1100 44.4159204
237 1200 57.9005269
238 1300 57.0619045
239 1400 66.1205699
240 1500 75.6204362
241 1600 74.8793775
242 1700 84.2891345
243 1800 83.3223266
244 1900 84.4328447
245 2000 92.2704318
246 2100 96.6705120
247 2200 90.5756019
248 2300 86.6116323
249 2400 97.2501922
250 2500 95.2831440
-------------------------------------------
Sincerely,
Philippe Grosjean
...........]<(({?<...............<?}))><...............................
) ) ) ) ) __ __
( ( ( ( ( |__) | _
) ) ) ) ) | hilippe |__)rosjean
( ( ( ( ( Marine Biol. Lab., ULB, Belgium
) ) ) ) ) __
( ( ( ( ( |\ /| |__)
) ) ) ) ) | \/ |ariculture & |__)iostatistics
( ( ( ( (
) ) ) ) ) e-mail: phgrosje at ulb.ac.be or phgrosjean at sciviews.org
( ( ( ( ( SciViews project coordinator (http://www.sciviews.org)
) ) ) ) ) tel: 00-32-2-650.29.70 (lab), 00-32-2-673.31.33 (home)
( ( ( ( (
) ) ) ) ) "I'm 100% confident that p is between 0 and 1"
( ( ( ( ( L. Gonick & W. Smith (1993)
) ) ) ) )
.......................................................................
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Possibly Parallel Threads
Wisdom of the Ancients
