1 问题的提出
问题:一批植物材料,设定了不同梯度的温度,每个处理重复3次,并测定了相对电导率,想计算各种基因型的半致死温度。
方法:采用文献《应用Logistic方程确定植物组织低温半致死温度的研究》的方法进行计算。
前言:很遗憾,R语言自带的nls函数(nls(ddl ~ SSlogis(temp, a, b, c))无法运行,主要缘故是数据问题。因此,根据查询的文献,从手动开始计算,到编程一次性完成。
2 文献所述求解方法
文献涉及的Logistic方程如下:
y
=
k
/
(
1
+
a
∗
e
(
−
b
∗
x
)
)
y=k/(1+a*e^{(-b*x)})
y=k/(1+a∗e(−b∗x))
问题改为涉及对参数
k
,
a
和
b
k,a和b
k,a和b的求解,半致死温度的求解公式如下:
L
t
50
=
l
n
a
/
b
Lt_{50}=ln{a}/b
Lt50=lna/b
关键步:参考文献需要对原始的电导率作如下转换:
y
′
=
(
k
−
y
)
/
y
=
a
∗
e
(
−
b
∗
x
)
y'=(k-y)/y=a*e^{(-b*x)}
y′=(k−y)/y=a∗e(−b∗x)
上述方程进行对数转换:
l
n
y
′
=
l
n
a
−
b
∗
x
ln{y'}=ln{a}-b*x
lny′=lna−b∗x
进一步
l
n
y
′
=
l
n
(
k
−
y
)
/
y
ln{y'}=ln{(k-y)/y}
lny′=ln(k−y)/y
式中,
y
y
y是原始的相对电导率(去掉%的值),
k
=
100
k=100
k=100。
通过
l
n
y
′
ln{y'}
lny′即可对
x
x
x做普通的线性回归,获得
l
n
a
和
b
ln{a}和b
lna和b的值。
最后要求解
k
k
k,采用等距离的3个
y
y
y值求解,方法如下:
k
=
(
y
2
2
∗
(
y
1
+
y
3
)
−
2
∗
y
1
∗
y
2
∗
y
3
)
/
(
y
2
2
−
y
1
∗
y
3
)
k=(y_{2}^2*(y_{1}+y_{3})-2*y_{1}*y_{2}*y_{3})/(y_{2}^2-y_{1}*y_{3})
k=(y22∗(y1+y3)−2∗y1∗y2∗y3)/(y22−y1∗y3)
通过上述的步骤,即可获得Logistic方程的各项参数值。之后就是利用R语言编程的事了,以及如何实现批量分析的目的。
3 R语言编程求解
3.1 数据读取和变换
对原始的电导率数据做所需的对数转换和均值计算:
df2<-readxl::read_excel('抗逆性.xls',sheet=1,skip=1)
df2$y<-log((100-df2$ddl)/df2$ddl)
df2a<-df2 %>% select(Clone,y,ddl,temp) %>%
group_by(Clone,temp) %>% summarise(dm=mean(ddl),ym=mean(y))
3.2 批量求解的函数
lgssf <- function(dat,mod=ym~temp,x=temp,y0=dm,tn=3:5,ts=3) {
# mod: line regression
# x: temperature, for plot
# y0: orginal leakage, for plot
# tn: used orginal leakage data, 3 points position
# ts: orginal leakage variable
x=deparse(substitute(x))
y0=deparse(substitute(y0))
m2<-lm(mod,data=dat)
tt<-summary(m2)
r2<-tt$r.squared
Fv<-tt$fstatistic
pv<-1-pf(Fv[1],df1=Fv[2],df2=Fv[3])
fit<-list()
fit$r2<-r2;fit$Fv<-Fv;fit$pv<-pv
ln.a<-coef(m2)[1];a<-exp(ln.a)
b<--coef(m2)[2]
Lt50<-ln.a/b
y1=dat[tn[1],ts];y2=dat[tn[2],ts];y3=dat[tn[3],ts]
k=(y2^2*(y1+y3)-2*y1*y2*y3)/(y2^2-y1*y3)
lgss.cof<-c(k,a,b,Lt50)
lgss.cof<-unlist(lgss.cof)
names(lgss.cof)<-c('k','a','b','Lt50')
lc<-lgss.cof<-round(lgss.cof,3)
equat<-paste0('y=',lc[1],'/','(1+',lc[2],'*exp^','(',lc[3],'*x))')
fig<-ggplot(dat, aes_string(x=x,y=y0))+geom_point()+
geom_smooth(method='auto',se=F)
res<-list(equat,lc[4],lgss.cof,fit)
names(res)<-c('equat','Lt50','coef','fit')
res$fig<-fig
return(res)
}
3.3 运行程序和查看结果
res2a <- plyr::dlply(df2a,"Clone",lgssf)
> names(res2a) # 各植物材料的基因型代码
[1] "1" "2" "3" "4" "5" "6" "7" "8" "9"
# 第一个基因型的结果
> res2a$`1`[c('equat','Lt50','fit')]
$`equat`
[1] "y=21.713/(1+1.148*exp^(-0.127*x))"
$Lt50
Lt50
-1.089
$fit
$fit$`r2`
[1] 0.8356973
$fit$Fv
value numdf dendf
20.34531 1.00000 4.00000
$fit$pv
value
0.01073532
图形查看res2a$`1`$fig
:
到此,所有求解和批量运算都结束了。