Hi everybody,
I'm using the method described here to make a linear regression: http://www.apsnet.org/education/advancedplantpath/topics/Rmodules/Doc1/05_Nonlinear_regression.html > ## Input the data that include the variables time, plant ID, and severity > time <- c(seq(0,10),seq(0,10),seq(0,10)) > plant <- c(rep(1,11),rep(2,11),rep(3,11)) > > ## Severity represents the number of > ## lesions on the leaf surface, standardized > ## as a proportion of the maximum > severity <- c( + 42,51,59,64,76,93,106,125,149,171,199, + 40,49,58,72,84,103,122,138,162,187,209, + 41,49,57,71,89,112,146,174,218,250,288)/288 > data1 <- data.frame( + cbind( + time, + plant, + severity + ) + ) > > ## Plot severity versus time > ## to see the relationship between## the two variables for each plant > plot( + data1$time, + data1$severity, + xlab="Time", + ylab="Severity", + type="n" + ) > text( + data1$time, + data1$severity, + data1$plant + ) > title(main="Graph of severity vs time") > > getInitial( + severity ~ SSlogis(time, alpha, xmid, scale), + data = data1 + ) alpha xmid scale 2.212468 12.506960 4.572391 > > > ## Using the initial parameters above, > ## fit the data with a logistic curve. > para0.st <- c( + alpha=2.212, + beta=12.507/4.572, # beta in our model is xmid/scale + gamma=1/4.572 # gamma (or r) is 1/scale + ) > > fit0 <- nls( + severity~alpha/(1+exp(beta-gamma*time)), + data1, + start=para0.st, + trace=T + ) 0.1621433 : 2.2120000 2.7355643 0.2187227 0.1621427 : 2.2124095 2.7352979 0.2187056 > > ## Plot to see how the model fits the data; plot the > ## logistic curve on a scatter plot > plot( + data1$time, + data1$severity, + type="n" + ) > > text( + data1$time, + data1$severity, + data1$plant + ) > > title(main="Graph of severity vs time") > > curve( + 2.21/(1+exp(2.74-0.22*x)), + from=time[1], + to=time[11], + add=TRUE + ) As you can see I have to do some work manually, such as setting the numbers to be used for calculation of alpha, beta and gamma. I wonder if you might have an idea how to automatize this? I suppose it should be possible to save the output from getInitial() and reach the elements via index or something, but how? I guess a similar approach could be used for the values of fit0? Or even better, if the variables alpha, beta and gamma could be used right away for instance in curve(), instead of adding the values manually. But just exchanging the values with the varables (alpha instead of 2.21 etc) doesn't seem to work. What is the reason for that? Any solution? A last, general but somewhat related question. If I set variables in a function such as para0.st <- c(alpha=2.212, ...), is it just stored locally, or can it be used globally, I mean, can I use the variable anywhere (for instance in curve()) or just in the function where it was created? I'm asking because I'm used to Java, where the life time of local variabels only extends to the closing braces, while global variables can be reached everywhere. The reason for automatization is that I'll have to repeat the procedure more than a hundred times, while making overview pair waise plots of my data, with both this logaritmic regression and several others (exponential, monomoelcular, logistic, Gompertz and Weibull). Wish you all the best, Joel _________________________________________________________________ Nya Windows 7 gör allt lite enklare. Hitta en dator som passar dig! http://windows.microsoft.com/shop [[alternative HTML version deleted]]
______________________________________________ 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.
