Sigmoid curve fitting for transpiration measurements from porometer at different water potentials (pressure):

por<-data.frame(list(structure(list(run = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L), .Label = c("1", "3", "4"), class = "factor"), press = c(15, 21, 24, 29.5, 15, 21, 24, 29.5, 15, 21, 24, 29.5), tr_rel = c(1, 0.459454191, 0.234697856, 0.135282651, 1, 0.853283066, 0.306314797, 0.186302231, 1, 0.42980063, 0.103882476, 0.086463799), tr = c(513, 235.7, 120.4, 69.4, 318.3, 271.6, 97.5, 59.3, 476.5, 204.8, 49.5, 41.2)), .Names = c("run", "press", "tr_rel", "tr"), row.names = c(NA, -12L), class = "data.frame"))) summary(mod1<-nls(tr ~ SSlogis( log(press), Asym, xmid, scal),data=por)) press_x <- seq(10, 40, length = 100) predict(mod1, data.frame(press = press_x)) with(por, plot(press,tr,xlim=c(10,35),ylim=c(0,500))) lines(press_x, predict(mod1, data.frame(press = press_x))) ## note that the variable in SSlogis is named "press" ## and the same name should appear on the left-hand side ## of "=" in the data frame ###http://finzi.psych.upenn.edu/R/Rhelp02/archive/42932.html ###calc. of conf. intervalls, by linear approximation se.fit <- sqrt(apply(attr(predict(mod1, data.frame(press = press_x)),"gradient"),1, function(x) sum(vcov(mod1)*outer(x,x)))) ###plot points and curve plus ci's matplot(press_x, predict(mod1, data.frame(press = press_x))+outer(se.fit,qnorm(c(.5, .025,.975))), type="l",lty=c(1,3,3),ylab="transp",xlab="waterpot") with(por, matpoints(press,tr,pch=1)) ###hhbio.wasser.tu-dresden.de/projects/modlim/doc/modlim.pdf ###R squared for non-linear Regression Rsquared <- 1 - var(residuals(mod1))/var(por$tr) ###add r^2 to plot: text(35,550,bquote(R^2==.(round(Rsquared, 3)))) asym<-coef(mod1)[1] y_tenth<-asym*0.1 abline(h=asym) abline(h=y_tenth,lty=5) text(10,65,"10% of max. transpiration",adj=0) ###finding y by x x<-15 y<-SSlogis(log(x), coef(mod1)[1], coef(mod1)[2] ,coef(mod1)[3]) y ###alternative way ###y = Asym/(1+exp((xmid-log(x))/scal)) y<-coef(mod1)[1]/(1+exp((coef(mod1)[2]-log(x))/coef(mod1)[3])) ###finding x by y Asym=coef(mod1)[1] xmid=coef(mod1)[2] scal=coef(mod1)[3] y_tenth<-asym*0.1 x_tenth<-exp((xmid-log(Asym/y_tenth-1)*scal));names(x_tenth)<-"x" x_tenth abline(v=x_tenth) text(10,25,paste("at waterpot.:", round(x_tenth,2)),adj=0)

Hi

ReplyDeletegood job.. Can you show me the sigmoidal equation. I'm going to apply this model to my data

thank you in advanced

Thank you very much. But I have two question:

Delete1) I just put my data and after "summary" appear this error "Error in nls(y ~ 1/(1 + exp((xmid - x)/scal)), data = xy, start = list(xmid = aux[1L], :

increasing factor 0.000488281 below 'minFactor' 0.000976562".....what is the meaning? I can't apply my data? Also with other kind of equation that are present in ?nls, it doesn't work. But if I exchange the data of x with y the equation works. What is the meaning?

2) If I want to find x from y in this script there is the equation "x_tenth=exp((xmid-log(Asym/y_tenth-1)*scal)). I don't understand the meaning of y_tenth. I mean if I put one value of waterpot (25 for example) I can find x with this equation? x=exp((xmid-log(Asym/25-1)*scal))

I'm sorry for all this questions, but I lost a lot of time without results

Thank you in advance

Its all there, see line 54 and 60!

DeleteHi

ReplyDeletei found very useful this script - which runs very well

i would like to apply it with an analysis involving SSfpl

here is the first part of my code

##

summary(mod1<-nls(DTest ~ SSfpl(ltc, a,b,c,d),data=dataMRr))

ltc_x <- seq(71, 224, length = 100)

predict(mod1, data.frame(ltc = ltc_x))

with(dataMRr, plot(ltc,DTest,xlim=c(70,230),ylim=c(0,40)), cex=0.2)

lines(ltc_x, predict(mod1, data.frame(ltc = ltc_x)), lwd=2)

se.fit <- sqrt(apply(attr(predict(mod1, data.frame(ltc = ltc_x)),"gradient"),1,

function(x) sum(vcov(mod1)*outer(x,x))))

##

which works well without any error message

but when i run

##

matplot(ltc_x, data.frame(ltc = ltc_x))+outer(se.fit,qnorm(c(.5, .025,.975)), type="l", lty=c(1,3,3))

##

i have the following message:

##

Error in outer(se.fit, qnorm(c(0.5, 0.025, 0.975)), type = "l", lty = c(1, :

using ... with FUN = "*" is an error

##

Any help would be greatly appreciated

This is old code - there might be some issues with package updates! Please post to R-Help or suchlike with your question and example code!

DeleteOK - many thanks Kay - i'll try on R-help

Deletecheers

Hi

ReplyDeleteI hope this tread is still active. I get the issue: "Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :

NA/NaN/Inf in 'x'"

if I only have 0's associated with a particular x value, anyway to get round this?

Division by zero!! See http://stackoverflow.com/questions/19544942/how-to-get-rid-of-error-with-multiple-regression-in-r

Delete