# Nonlinear Regression / Binding Curve
# simulate noisy data for myoglobin KD = 0.26 [L] in kPa
# includes example of a loop
# data generated with random noise added
# Start
KD <- 0.26 # initializing variables
L <- 0.0
y <- 0.0
for (i in 1:15) { y[i] <- (0.10*i/(KD+0.10*i)) + rnorm(1, sd = 0.03)
L[i] <- 0.10*i }
L
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
y
## [1] 0.2788516 0.4285463 0.5512855 0.5790450 0.6619387 0.6880949 0.7105009
## [8] 0.7026543 0.7530362 0.7943787 0.8025836 0.8465140 0.7950814 0.8661941
## [15] 0.8771812
# try to do this with vectorized code!
# this is error
i <- seq(.1,1.5,0.1)
i
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
y <- i/(KD+i) + rnorm(1, sd = 0.03)
y
## [1] 0.3013756 0.4583804 0.5593121 0.6296584 0.6814926 0.7212723 0.7527645
## [8] 0.7783148 0.7994599 0.8172486 0.8324214 0.8455156 0.8569312 0.8669713
## [15] 0.8758706
tryfit <- nls(y ~ L/(KD+L),
start = c(KD = 0.5))
summary(tryfit)
##
## Formula: y ~ L/(KD + L)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## KD 0.229997 0.001794 128.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005541 on 14 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 1.141e-06
plot(L,y, main = "Myoglobin Binding Curve", xlab = "pO2 Kpa", ylab = "Fractional Binding")
lines(L,predict(tryfit), col = "blue") # add a line fitted to points
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.