R语言青霉素稳定性研究编程简略代码

刘先聚

tempr<-c(37,43,54,60)
k<-c(NA,NA,NA,NA)
datals<-data.frame(tempr,k)

temp<- c(rep(tempr[1],6),rep(tempr[2],6),rep(tempr[3],6),rep(tempr[4],6))   #温度
time<-c(0,2,5,10,15,20,                     #取样时间
        0,1,2,4,8,12,
        0,0.5,1,2,3,4,
        0,0.25,0.5,1,1.5,2)
concen<-c(100,96.03,88.18,80.574,72.32,64.92,      #浓度
          100,96.025,92.256,85.08,71.433,62.305,
          100,94.035,87.962,79.31,69.016,61.392,
          100,94.94,89.067,82.36,73.308,64.13)
datax<-data.frame(temp,time,concen)
for(i in 1:4){
newdata<-datax[datax$temp==tempr,]
with(newdata,{
y<-concen
x<-time
nlmod<-nls(y~A*exp(-B*x))
datals[i,2]<<-coefficients(nlmod)
})
}
x<-1/(datals$tempr+273.15)
y<-log(datals$k)
datalz<-data.frame(y,x)
datalz
lmod<-lm(y~x,data=datalz)
summary(lmod)
coefficients(lmod)
u<-coefficients(lmod)
k25<-exp(u[1]+u[2]/(25+273.15))
t0.9<- -log(0.9)/k25
t0.9datax
   temp  time  concen
1    37  0.00 100.000
2    37  2.00  96.030
3    37  5.00  88.180
4    37 10.00  80.574
5    37 15.00  72.320
6    37 20.00  64.920
7    43  0.00 100.000
8    43  1.00  96.025
9    43  2.00  92.256
10   43  4.00  85.080
11   43  8.00  71.433
12   43 12.00  62.305
13   54  0.00 100.000
14   54  0.50  94.035
15   54  1.00  87.962
16   54  2.00  79.310
17   54  3.00  69.016
18   54  4.00  61.392
19   60  0.00 100.000
20   60  0.25  94.940
21   60  0.50  89.067
22   60  1.00  82.360
23   60  1.50  73.308
24   60  2.00  64.130

25℃时 t.0.9= 18.39h

syhorchid

编程啊，不懂哦。

刘先聚

syhorchid 发表于 2019-3-13 23:56
编程啊，不懂哦。

下载一个RStudio软件
然后将代码部分贴进去就可以计算出t0.9了本质的部分还是阿累尼乌斯方程的应用
如果有其它数据只要修改代码部分的温度 时间 浓度的数据即可。
不过一般很少用而已

刘先聚

# 创建温度向量
tempr <- c(37, 43, 54, 60)
# 创建 NA 值的向量，长度与 tempr 相同
k <- c(NA, NA, NA, NA)
# 创建数据框 datals，包含 tempr 和 k
datals <- data.frame(tempr, k)

# 重复 tempr 中的值，使得每个温度对应六个时间点
temp <- c(rep(tempr[1], 6), rep(tempr[2], 6), rep(tempr[3], 6), rep(tempr[4], 6))
# 创建时间向量 time，对应不同温度下的取样时间
time <- c(0, 2, 5, 10, 15, 20,  # 37°C 下的时间点
          0, 1, 2, 4, 8, 12,     # 43°C 下的时间点
          0, 0.5, 1, 2, 3, 4,    # 54°C 下的时间点
          0, 0.25, 0.5, 1, 1.5, 2)# 60°C 下的时间点
# 创建浓度向量 concen，对应不同温度下的浓度值
concen <- c(100, 96.03, 88.18, 80.574, 72.32, 64.92,  # 37°C 下的浓度
            100, 96.025, 92.256, 85.08, 71.433, 62.305, # 43°C 下的浓度
            100, 94.035, 87.962, 79.31, 69.016, 61.392, # 54°C 下的浓度
            100, 94.94, 89.067, 82.36, 73.308, 64.13)   # 60°C 下的浓度
# 创建数据框 datax，包含温度、时间和浓度
datax <- data.frame(temp, time, concen)

# 循环处理每个温度下的数据
for (i in 1:4) {
  # 选择当前温度下的数据
  newdata <- datax[datax$temp == tempr[i], ]
  # 使用 with 函数简化 newdata 的调用
  with(newdata, {
    # 指定浓度和时间作为变量
    y <- concen
    x <- time
    # 使用非线性最小二乘法拟合数据
    nlmod <- nls(y ~ A * exp(-B * x))
    # 从拟合模型中提取系数 B
    datals[i, 2] <<- coefficients(nlmod)["B"]
  })
}

# 输出 datals 数据框
print(datals)

# 计算绝对温度的倒数
x <- 1 / (datals$tempr + 273.15)
# 计算 k 的自然对数
y <- log(datals$k)
# 创建新数据框 datalz
datalz <- data.frame(y, x)

# 使用线性模型拟合 y ~ x
lmod <- lm(y ~ x, data = datalz)

# 输出线性模型的摘要信息
summary(lmod)

# 输出线性模型的系数
print(coefficients(lmod))

# 提取线性模型的系数
u <- coefficients(lmod)

# 计算 25°C 下的 k 值
k25 <- exp(u[1] + u[2] / (25 + 273.15))

# 计算浓度降至 90% 所需的时间 t0.9
t0.9 <- -log(0.9) / k25

# 输出 t0.9 的值
print(t0.9)

# 输出 datax 数据框
print(datax)

刘先聚

> print(datals)
  tempr         k
1   37 0.04501756
2   43 0.05765624
3   54 0.12691752
4   60 0.18771560

Call:
lm(formula = y ~ x, data = datalz)

Residuals:
        1         2         3         4
0.0017905 0.0005701 0.0014901 0.0011492

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  4.15552    0.01101  377.42 2.07e-05 ***
x           -4.59460    0.13040  -35.24 2.74e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001574 on 2 degrees of freedom
Multiple R-squared:  0.9984,        Adjusted R-squared:  0.9979
F-statistic:  1240 on 1 and 2 DF,  p-value: 2.738e-04

> print(coefficients(lmod))
(Intercept)           x
4.15552025 -4.59459967

> print(t0.9)
[1] 18.38992

> print(datax)
   temp  time  concen
1    37  0.00 100.000
2    37  2.00  96.030
3    37  5.00  88.180
4    37 10.00  80.574
5    37 15.00  72.320
6    37 20.00  64.920
7    43  0.00 100.000
8    43  1.00  96.025
9    43  2.00  92.256
10   43  4.00  85.080
11   43  8.00  71.433
12   43 12.00  62.305
13   54  0.00 100.000
14   54  0.50  94.035
15   54  1.00  87.962
16   54  2.00  79.310
17   54  3.00  69.016
18   54  4.00  61.392
19   60  0.00 100.000
20   60  0.25  94.940
21   60  0.50  89.067
22   60  1.00  82.360
23   60  1.50  73.308
24   60  2.00  64.130