Re: [R] nls: Fitting two models at once?

From: Heather Turner <Heather.Turner_at_warwick.ac.uk>
Date: Fri, 07 Nov 2008 10:46:24 +0000

Dear Martin,

You can use the same idea of concatenating the data in R. The following reproduces your example:

spec <- c(asfe, dias)
n1 <- length(asfe)
n2 <- length(dias)
CYT <- c(28 * CYTF, 24.6 * 2 * CYTB6)
B559HP <- c(B559HP, numeric(n2))
B559LP <- c(numeric(n1), B559LP)
C550 <- c(numeric(n1), C550)

spec.fit <-
nls(

	spec ~ cyt.v*CYT + hp.v*B559HP + lp.v*B559LP + c550.v*C550,
	start = list(cyt.v = 0.00005, hp.v = 0.000003, lp.v = 1,
	             c550.v = 1) # arbitrary

);

or more efficiently, using lm

spec.fit <- lm(spec ~ 0 + CYT + B559HP + B559LP + C550)

# draw stuff
plot(

	1, 2,
	type="n",
	xlim = c(540, 575),
	ylim=c(-0.002, 0.008),

);

# first spectrum and fit
lines(wl, asfe, type="b", pch=19); # solid circles lines(wl, fitted(spec.fit)[1:n1], col = "red");

# second spectrum and fit
lines(wl, dias, type="b");
lines(wl, fitted(spec.fit)[(n1 + 1):(n1 + n2)], col = "blue");

Best wishes,

Heather

Dr H Turner
Senior Research Fellow
Dept. of Statistics
The University of Warwick
Coventry
CV4 7AL

Tel: 024 76575870
Fax: 024 76524532
Url: www.warwick.ac.uk/go/heatherturner


Martin Ballaschk wrote:
> Hello,
>
> I'm still a newbie user and struggling to automate some analyses from
> SigmaPlot using R. R is a great help for me so far!
>
> But the following problem makes me go nuts.
>
> I have two spectra, both have to be fitted to reference data. Problem: the
> both spectra are connected in some way: the stoichiometry of coefficients
> "cytf.v"/"cytb.v" is 1/2.
> {{In the SigmaPlot workflow one has to copy the two spectra into one column
> beneath each other and the two spectra are somehow treated as one curve -
> like in
http://home.arcor.de/ballaschk/cytbf-help/sigmaplot%20formula.png}}
>
> Can anybody help? :(
>
> I tried to condense everything to the "minimum" R script below to give an
> impression of what I'm talking about.
>
> Martin
>
>
>
> ########################################################################
> # "Minimal" R script reading remote data for convenience
>
> ### READ IN DATA
> # first spectrum
> asfe <- read.table("http://home.arcor.de/ballaschk/cytbf-help/asfe.csv")[, 1];
> # second spectrum
> dias <- read.table("http://home.arcor.de/ballaschk/cytbf-help/dias.csv")[, 1];
>
> # reference data for fit, wavelength = wl
> ref <- read.table("http://home.arcor.de/ballaschk/cytbf-help/reference.csv",
> sep="\t", dec=".", header=T);
> attach(ref);
>
> ### FITTING, problem: 2*cytf.v == cytb.v
>
> # fit first spectrum to two reference spectra
> asfe.fit <-
> nls(
> asfe ~ (cytf.v * 28) * CYTF + hp.v * B559HP,
> start = list(cytf.v = 0.00005, hp.v = 0.000003) # arbitrary
> );
>
> # fit second spectrum to three reference spectra
> dias.fit <-
> nls(
> dias ~ (cytb.v * 24.6 * 2) * CYTB6 + lp.v * B559LP + c550.v * C550,
> start = list(cytb.v = 1, lp.v = 1, c550.v = 1) # arbitrary
> );
>
>
> # draw stuff
> plot(
> 1, 2,
> type="n",
> xlim = c(540, 575),
> ylim=c(-0.002, 0.008),
> );
>
> # first spectrum and fit
> lines(wl, asfe, type="b", pch=19); # solid circles
> lines(wl, fitted(asfe.fit), col = "red");
>
> # second spectrum and fit
> lines(wl, dias, type="b");
> lines(wl, fitted(dias.fit), col = "blue");
>
> ______________________________________________
> R-help_at_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.



R-help_at_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. Received on Fri 07 Nov 2008 - 10:52:17 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 07 Nov 2008 - 22:30:24 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive