批量求解基于电导率的植物半致死温度

批量求解基于电导率的植物半致死温度

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 variablex=deparse(substitute(x))y0=deparse(substitute(y0))m2<-lm(mod,data=dat)tt<-summary(m2)r2<-tt$r.squaredFv<-tt$fstatisticpv<-1-pf(Fv[1],df1=Fv[2],df2=Fv[3])fit<-list()fit$r2<-r2;fit$Fv<-Fv;fit$pv<-pvln.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*f<-c(k,a,b,f<-f)f)<-c('k','a','b','Lt50')lc<-f<-f,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],f,fit)names(res)<-c('equat','Lt50','coef','fit')res$fig<-figreturn(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))"$Lt50Lt50 
-1.089 $fit
$fit$`r2`
[1] 0.8356973
$fit$Fvvalue    numdf    dendf 
20.34531  1.00000  4.00000
$fit$pvvalue 
0.01073532

图形查看res2a$`1`$fig：

到此，所有求解和批量运算都结束了。