hello,
i managed to get CI's for my curve - but now I need the intervall for a
certain y point (y_tenth) of the curve..
can anyone help me with this?
#####data:
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)))
###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 with 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))
###the desired y is the tenth of the y asymptote
asym<-coef(mod1)[1]
tenth<-asym*0.1
abline(h=asym)
abline(h=tenth,lty=5)
text(10,65,"10% of max. transpiration",adj=0)
###finding corresponding x by y with inverse function:
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)
--
View this message in context:
http://n4.nabble.com/getting-CI-s-for-certain-y-of-nls-fitted-curve-tp1695025p1695025.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
