0

This is a similar question to that posted in Regression (logistic) in R: Finding x value (predictor) for a particular y value (outcome). I am trying to find the x value for a known y value (in this case 0.000001) obtained from fitting a log normal curve fitted to sapling densities at distances from parent trees using a genetic algorithm. This algorithm gives me the a and b parameters of the best-fit log normal curve.

I have obtained the value of x for y=0.00001 for other curves, such as negative exponential, by using uniroot using this code (which works well for these curves):

##calculate x value at y=0.000001 (predicted near-maximum recruitment distance)
aparam=a
bparam=b
testfn <- function (y, aparam, bparam) {
## find value of x that satisfies y = a + bx
    fn <- function(x) (a * exp(-b * x)) - y
    uniroot(fn, lower=0, upper= 100000000)$root
}
testfn(0.000001)

Unfortunately, the same code using a log normal formula does not work. I have tried to use uniroot by setting the lower boundary above zero. But get an error code:

Error in uniroot(fn, lower = 1e-16, upper = 1e+18) : 
  f() values at end points not of opposite sign

My code and data (given below the code) is:

file="TR maire 1mbin.txt"
xydata <- read.table(file,header=TRUE,col.names=c('x','y'))



####assign best parameter values
a = 1.35577
b = 0.8941521



#####Plot model against data
par(mar=c(5,5,2,2))

xvals=seq(1,max(xydata$x),1)
plot(jitter(xydata$x), jitter(xydata$y),pch=1,xlab="distance from NCA (m)",
ylab=quote(recruit ~ density ~ (individuals ~ m^{2~~~ -1})))
col2="light grey"


plotmodel <-  a* exp(-(b) * xvals)
lines(xvals,plotmodel,col=col2)




####ATTEMPT 1
##calculate x value at y=0.000001 (predicted near-maximum recruitment distance)
aparam=a
bparam=b
testfn <- function (y, aparam, bparam) {
    fn <- function(x) ((exp(-(((log(x/b)) * (log(x/b)))/(2*a*a))))/(a * x * sqrt(2*pi))) - y
    uniroot(fn, lower=0.0000000000000001, upper= 1000000000000000000)$root
}
testfn(0.000001)

data is:

 xydata

1    1 0.318309886
2    2 0.106103295
3    2 0.106103295
4    2 0.106103295
5    3 0.063661977
6    4 0.045472841
7    5 0.035367765
8    5 0.035367765
9    7 0.048970752
10   8 0.021220659
11   8 0.021220659
12   8 0.042441318
13   9 0.018724111
14  10 0.016753152
15  10 0.016753152
16  12 0.013839560
17  13 0.025464791
18  16 0.010268061
19  17 0.009645754
20  24 0.013545102
21  25 0.032480601
22  26 0.043689592
23  27 0.006005847
24  28 0.011574905
25  31 0.062618338
26  32 0.005052538
27  42 0.003835059
28  42 0.003835059
29  44 0.003658734
30  46 0.003497911
31  48 0.006701261
32  50 0.003215251
33  50 0.006430503
34  51 0.006303166
35  58 0.002767912
36  79 0.002027452
37 129 0.003715680
38 131 0.001219578
39 132 0.001210304
40 133 0.001201169
41 144 0.001109094
42 181 0.000881745
43 279 0.001142944
44 326 0.000488955

Or is there another way of approaching this? I'm an ecologist and sometimes R just does not make sense!

IRTFM
  • 258,963
  • 21
  • 364
  • 487
Mike Thorsen
  • 1
  • Possible duplicate of [Uniroot solution in R](https://stackoverflow.com/questions/38961221/uniroot-solution-in-r) – user10089632 Aug 28 '17 at 23:00
  • Yes, I saw that post too. I couldn't plot my testfn to see its shape and didn't see what the solution to the problem is. – Mike Thorsen Aug 28 '17 at 23:26
  • 1
    learning demands some effort, at your place plotting my function is the first thing I do. Plus, if you have attentively read that post, you would understand that this is the expected behaviour – user10089632 Aug 28 '17 at 23:31
  • Fair enough. I have been working on this for over a day now. I can't even plot test fn: – Mike Thorsen Aug 28 '17 at 23:40
  • Edit the above as timed out: Fair enough. I have been working on this for over a day now. I can't even plot test fn: lines(xvals,testfn) as I get the error code "Error in xy.coords(x, y) : 'x' and 'y' lengths differ". Given the formula I'm not sure what x and y even are. I just get more and more confused. – Mike Thorsen Aug 28 '17 at 23:47

1 Answers1

0

Seems like there were some errors in my r code, but the main problem is that my lower limit was too low and the Log Normal curve does not extend to that value (my interpretation). The solution that works for me is:

    ### define the formula parameter values

        a = 1.35577
        b = 0.8941521


        ### define your formula (in this instance a log normal) in the {}
        fn <- function(x,a,b,y) { ((exp(-(((log(x/b)) * (log(x/b)))/(2*a*a))))/(a * x * sqrt(2*pi))) - y}


###then use uniroot()$root calling the known parameter values and defining the value of y that is of interest (in this case 0.000001)

        uniroot(fn,c(1,200000),a=a,b=b,y=0.000001)$root
Mike Thorsen
  • 1