我花了两个晚上很专心的阅读,总算把一篇文章(Watkins P. and Venables B. 2006. Non-linear regression for optimizing the separation of carboxylic acids. R News 6(3): 2-7. )的每一个R命令都弄得差不多明白了(自认为,还需要大家确认!)。最后发现确实值得深入学习一下。除了进一步理解了对nls的self-starter function的编写之外,还学到了一些非常好用,而且又能够在今后经常用到的functions, 比如 eval(), with(), matplot(), matpoints(), transform(),以及对nls的update().
整个阅读的过程伴随着把文中的命令输入R跑一下,然后看帮助,然后把一些小的部分拆解开跑着看等方法,现在把命令贴在下面,详细的学习备忘,改日贴上。
# input data <br />
don<-data.frame(ph=c(3.79,3.79,4.14,4.38,4.57,4.74,4.74,4.92,5.11,5.35,5.67,5.67), <br />
ba=c(34.21,34.27,25.85,20.46,15.61,12.42,11.42,9.64,7.3,5.15,3.18,3.18), <br />
oaba=c(15.06,14.64,14.24,13.33,12.61,11.33,10.55,10.15,9.12,6.36,3.92,3.92), <br />
paba=c(8.85,8.33,8.00,7.58,6.82,5.76,5.76,5.09,4.15,2.88,1.6,1.58), <br />
hoba=c(14.3,14.52,12.3,10.76,8.91,7.24,7.06,5.94,4.52,3.09,1.68,1.62)) <br />
<br />
# transform pH to [+H] <br />
tmp<-transform(don,H=10^(-ph)) <br />
# 非常怪异的方法,我费了不少时间去弄明白文章每句话的意思,详见后。 <br />
<br />
##### so called "partially linear" models ###### <br />
ba.nls<-nls(ba~cbind(1,H/ka)/(1+H/ka),data=tmp,algorithm="plinear",start=c(ka=0.0001),trace=TRUE) <br />
<br />
## the trace display <br />
#13.24040 : 0.000100 3.534884 54.813991 <br />
#7.178175 : 4.295106e-05 5.944938e-01 4.197216e+01 <br />
#1.439102 : 5.602806e-05 1.640794e+00 4.525283e+01 <br />
#1.274617 : 5.908893e-05 1.836970e+00 4.597673e+01 <br />
#1.274311 : 5.923005e-05 1.845664e+00 4.600979e+01 <br />
#1.274311 : 5.923158e-05 1.845757e+00 4.601015e+01 <br />
#估计得到的参数 <br />
coef(ba.nls) <br />
#看一下拟合效果 <br />
plot(don$ph,don$ba,ylim=c(0,40)) <br />
lines(don$ph,fitted(ba.nls),col=2) <br />
<br />
##############Self-starting function############## <br />
##作者说 The initial value routine has to use a somewhat obscure convention way. <br />
## This convention initially appears to be obscure, but it becomes less so with time! <br />
##可见我觉得晦涩有情可原! <br />
<br />
# Initial value routine <br />
SSba.init<-function(mCall, data,LHS){ <br />
# <br />
#k<-(k1+k0*H/ka)/(1+H/ka); H=10^(-ph) <br />
# <br />
<br />
H<-10^(-eval(mCall[["ph"]],data)) <br />
k<-eval(LHS,data) <br />
ka<-as.vector(coef(lsfit(cbind(H,-k),H*k,int=TRUE))[3]) <br />
b<-coef(nls(k~cbind(1,H/ka)/(1+H/ka),data=data.frame(k=k,H=H), <br />
algorithm="plinear", <br />
start=c(ka=ka))) <br />
names(b)<-mCall[c("ka","k1","k0")] <br />
b <br />
} <br />
<br />
#To understand the calcalation of ka <br />
#?lsfit() # least square estimate <br />
don<-transform(don,H=10^-(-ph)) <br />
xx<-with(don,cbind(H,-ba)) <br />
yy<-with(don,H*ba) <br />
ls.1<-lsfit(xx,yy,int=TRUE) # int==intercept <br />
coef(ls.1) # the order is "intercept", "H","K" <br />
##Then I understand why it is "coef(ls.1)[3]" <br />
<br />
<br />
<br />
<br />
## self-starting model <br />
<br />
SSba<-selfStart(~(k1+k0*10^(-ph)/ka)/(1+10^(-ph)/ka), <br />
initial=SSba.init, <br />
parameters=c("ka","k1","k0"), <br />
template=function(ph,k1,k0,ka){} <br />
) <br />
<br />
#run the model <br />
ba.ss<-nls(ba~SSba(ph,k1,k0,ka),data=don,trace=TRUE) <br />
<br />
# check the fitting (the same as "partially linear") <br />
lines(don$ph,fitted(ba.ss),col=4) <br />
<br />
# note the final "point"! <br />
oaba.ss<-update(ba.ss,oaba~.) <br />
paba.ss<-update(ba.ss,paba~.) <br />
hoba.ss<-update(ba.ss,hoba~.) <br />
<br />
<br />
points(don$ph,don$oaba,pch=2) <br />
lines(don$ph,fitted(oaba.ss),col=2,lty=2) <br />
<br />
## display the fitting <br />
form<-Quote((k1+k0*10^(-ph)/ka)/(1+10^(-ph)/ka)) <br />
<br />
at<-function(x,m) c(x,as.list(coef(m))) <br />
phdata<-list(ph=seq(3.5,6,len=250)) <br />
phdata<-transform(phdata,ba=eval(form,at(phdata,ba.ss)), <br />
ba=eval(form,at(phdata,ba.ss)), <br />
oaba=eval(form,at(phdata,oaba.ss)), <br />
paba=eval(form,at(phdata,paba.ss)), <br />
hoba=eval(form,at(phdata,hoba.ss)) <br />
) <br />
<br />
# plot the fitted lines <br />
with(phdata,matplot(ph,cbind(ba,oaba,paba,hoba),col=1:4,type="l",lty=1:4)) <br />
## points <br />
with(don,matpoints(ph,cbind(ba,oaba,paba,hoba),col=1:4,pch=1:4,cex=0.8)) <br />
ps:
1)多谢益辉帮我传了图片!
2)很高兴看到熟人陈沉, 为什么让俺注意休息?难道是看我发贴时间较晚? 估计是时差的缘故 [s:13]
3)突然意识到想说感谢“谢益辉”不容易[s:11]