Disclaimer
I am new to this site, relatively new to R (two weeks of learning), have just a really basic knowledge in statistics so sorry if I'm doing a dumb mistake there or asking bad question or something.
I also don't know how to nicely embed my dataset into post (I searched meta without finding anything about this) so I shared a link from Google Drive (if you know a better way, please let me know and I can change it or do it yourself if you have the right to edit).
Explanation
In physical Laboratory course I am attending I have done an experiment to determine attenuation coefficient of an optical medium using optical filters of different thicknesses (with this spectrometer http://www.vernier.com/products/sensors/spectrometers/visible-range/v-spec/).
It generated data which you can download here to reproduce my results: https://docs.google.com/uc?authuser=0&id=0B5x0REqWCRBuaHlXc2JWZ19sWE0&export=download
My goal is to determine $\kappa(\lambda)$, where $\kappa$ is the attenuation coefficient and $\lambda$ wavelenght, using Beer–Lambert law which states ($l$ is the thickness): $$\theta=x_0 \cdot \exp(-\kappa \cdot l)$$
I have 5 "curves" for different thicknesses of optical filters (from 1 to 5 mm) and I fit the exponential model for each wavelenght $\lambda$ which gives me the $\kappa$ for that wavelenght.
I tried this:
Linearized model
by taking logarithm of the equation above and some algebraic manipulation we get: $$ \ln \left( \frac{x_0}{\theta} \right) = \kappa \cdot l$$
represented by R command lm( I(log(x0/temp.theta)) ~ l + 0 )
Non-linear model
I take the Beer–Lambert law (first equation) and model it by nls(temp.theta ~ I(x0 * exp(-k*l)) + 0 , start= list(k = k.l)), where is obtained from the linearized model k.l = coef(lm(...)).
Code for reproduction
I set x0 = 100 since $\theta$ is measured in percents.
x0 = 100 #x=0 intercept of the exp function
data = read.delim2("export2.csv")
l = data$l
l = l[-which(is.na(l))]
data$l = NULL
data$kappa.l = NA
data$sd.l = NA
data$kappa.nl = NA
data$sd.nl = NA
for(i in 1:nrow(data)){
temp.theta = as.numeric(data[i,-which(names(data) %in% c("lambda","kappa.l","kappa.nl","sd.l","sd.nl"))])
temp.lm = lm( I(log(x0/temp.theta)) ~ l + 0 )
k.l = coef(temp.lm)
data[i , "kappa.l"] = k.l
data[i , "sd.l"] = coef(summary(temp.lm))[ ,"Std. Error"]
temp.nls = nls(temp.theta ~ I(x0 * exp(-k*l)) + 0, start = list(k = k.l))
data[i , "kappa.nl"] = coef(temp.nls)
data[i , "sd.nl"] = coef(summary(temp.nls))[ ,"Std. Error"]
}
Questions
When I visualize the result I get the graph below. These questions emerge:
- Why are the fits so much different around the 400 nm wavelenght?
- Which one fits the data better?
- Or are the data so far from exponential behavior there so I have to cut them where they meet?
- Which of these models is statistically (more) correct?
- Can I make the fits better by adding weights according to the measurement uncertainty (which I believe to be 5 % as stated on the Vernier web)?
the fit with higher values is the non-linear one
