1

I have some data generated using the following lines of code,

x <- c(1:10)
y <- x^3
z <- y-20
s <- z/3
t <- s*6
q <- s*y
x1 <- cbind(x,y,z,s,t,q)
x1 <- data.frame(x1)

I would like to plot x versus y,s, and t so I melt the data frame x1 first,

 library(reshape2)

 xm <- melt(x1, id=names(x1)[1], measure=names(x1)[c(2, 4, 5)], variable = "cols"`)

Then I plot them along with their linear fits using the following code,

library(ggplot2)
plt <- ggplot(xm, aes(x = x, y = value, color = cols)) +
  geom_point(size = 3) +
  labs(x = "x", y = "y") + 
  geom_smooth(method = "lm", se = FALSE)
plt

The plot which is generated is shown below,

enter image description here

Now I would liked to interpolate the x-intercept of the linear fit. The point in the plot where y axis value is 0.

The following lines of code as shown here, extracts the slope and y-intercept.

fits <- by(xm[-2], xm$cols, function(i) coef(lm(value ~ x, i)))
data.frame(cols = names(fits), do.call(rbind, fits))

Is there any way how I can extract the x-intercept other than manually calculating from the slope and y-intercept?

Thanks for the help!

Community
  • 1
  • 1
Amm
  • 1,749
  • 4
  • 17
  • 27

3 Answers3

3

You could do inverse prediction as implemented in package chemCal for calibrations if you don't want to calculate this yourself:

library(chemCal)
res <- by(xm[-2], xm$cols, function(i) inverse.predict(lm(value ~ x, i), 0)$Prediction)
res[1:3]
#xm$cols
#y        s        t 
#2.629981 2.819734 2.819734

Edit:

Maybe you prefer this:

library(plyr)
res <- ddply(xm, .(cols), 
  function(i) data.frame(xinter=inverse.predict(lm(value ~ x, i), 0)$Prediction))
#   cols   xinter
# 1    y 2.629981
# 2    s 2.819734
# 3    t 2.819734
Roland
  • 127,288
  • 10
  • 191
  • 288
  • Thanks for the suggestion. I would like to get the result as a single column, so I use `t(res[1:3])` which does not work. Here in this case there are only 3 fits so we give the index 1:3 but what if there are _n_ number of fits, where one does not know the value of _n_. – Amm Apr 16 '14 at 08:34
1

I don't think you can avoid computing the linear equation, though of course you don't have to do it by hand (unless you want to). For example:

by(xm[-2], xm$cols, function(i) { 
fit <- lm(value~x, i); print(fit); solve(coef(fit)[-1], -coef(fit)[1] )}
)

Call:
lm(formula = value ~ x, data = i)

Coefficients:
(Intercept)            x  
     -277.2        105.4  


Call:
lm(formula = value ~ x, data = i)

Coefficients:
(Intercept)            x  
     -99.07        35.13  


Call:
lm(formula = value ~ x, data = i)

Coefficients:
(Intercept)            x  
     -594.4        210.8  

xm$cols: y
[1] 2.629981
----------------------------------------------------------------------------------------------------------------- 
xm$cols: s
[1] 2.819734
----------------------------------------------------------------------------------------------------------------- 
xm$cols: t
[1] 2.819734

What was solved is basically -277.2 + 105.4*x = 0 for x -> 105.4*x = 277.2 (the solve-function call) -> x = 2.629981. Seems your lines 's' and 't' intersect the y=0 axis at the same spot. If I understood correctly, your problem isn't extrapolation since your x-range covers the intercept but instead interpolation.

Ps. I think your code was missing: require("reshape")

EDIT:

result <- c(by(xm[-2], xm$cols, function(i) { fit <- lm(value~x, i); print(fit); solve(coef(fit)[-1], -coef(fit)[1] )} )); print(result) 
> print(result)
       y        s        t 
2.629981 2.819734 2.819734
Teemu Daniel Laajala
  • 2,316
  • 1
  • 26
  • 37
  • Thanks for your suggestion. I have also found another method to solve this one now. I edited my code to include reshape package now. – Amm Apr 10 '14 at 12:33
  • is there a way to pass the values of the x-intercept to a data frame in the result list which was generated. – Amm Apr 10 '14 at 12:55
  • I'm not entirely sure what you mean... do you mean how to access the raw values returned by the solver? You could do e.g. : result <- c(by(xm[-2], xm$cols, function(i) { fit <- lm(value~x, i); print(fit); solve(coef(fit)[-1], -coef(fit)[1] )} )); print(result) – Teemu Daniel Laajala Apr 14 '14 at 23:34
  • I would like to get the result in a single column in the format which is indicated in my answer below. – Amm Apr 15 '14 at 06:38
  • Ok, I edited my answer to show this. With my answer it gives a row vector, but you can do "t(result)" to obtain the same column vector that you have in yours. – Teemu Daniel Laajala Apr 15 '14 at 12:23
  • I was thinking of transpose as well. Thanks for your suggestion. – Amm Apr 15 '14 at 12:42
0

I found a way to calculate the x-intercept, first create a data frame with the y-intercept and slope values,

par <- data.frame(cols = names(fits), do.call(rbind, fits))

Then rename column header names to accurately denote the values,

colnames(par)[2] <- "y_intercept"
colnames(par)[3] <- "slope"
# Calculate the x-intercept by using the formula -(y_intercept)/slope
x_incpt <- -par[2]/par[3]
colnames(x_incpt) <- "x_intercept"

Which gives the following result,

  x_intercept
y    2.629981
s    2.819734
t    2.819734
Amm
  • 1,749
  • 4
  • 17
  • 27