机器算法验证 - 如何在 R 中使用数学模型进行数据分析 - 吾爱随笔录

如何在 R 中使用数学模型进行数据分析

机器算法验证 r 回归生态

2022-04-10 07:50:50

我希望使用 Firbank & Watkinson (1985) J. App 开发的数学模型。生态。22:503-517 用于分析混合种植的植物之间的竞争。

模型如下：

W_{A} = W_{m A} {(1 + a_{A} (N_{A} + α N_{B}))}^{- b_{A}}

$W_{A}=W_{mA}\left(1 + a_{A}\left(N_{A}+\alpha N_{B}\right)\right)^{-b_{A}}$

其中是实验中种植每株植物的平均产量是物种的分离植物的平均产量，所需的表面积,是物种的种植密度，是竞争系数，是“资源利用效率”参数。 $W_{A}$ $A$ $W_{mA}$ $A$ $a_{A}$ $W_{mA}$ $N_{A}$ $A$ $\alpha$ $-b_{A}$

据我了解，该模型是回归模型。我有物种和的密度数据以及（和）以及响应变量。我不确定如何使用 R 来估计剩余值，其中最重要的是竞争系数。如果我需要提供更多信息，请告诉我。 $A$ $B$ $N_{A}$ $N_{B}$ $W_{A}$ $\alpha$

2个回答

您可以使用非线性回归将此方程拟合到您的数据中。

我会试试看nls。使用的关键方面nls是提供合理的起始值。示例代码可能类似于

nls.mod <- nls(Wa~Wma*(1+a*(Na + alpha*N))^(-b), data = dataset, start = list(a = 1, b = 1, alpha = 1)). – COOLSerdash

非线性回归也是我的第一个想法——但要注意，非线性最小二乘法默认假设方差不变；一旦您正确考虑了误差项，可能是等式的修改版本（例如可能是对数刻度）可能是对关系的更好描述。如果您没有理论作为选择误差项的一种方式，您可以查看数据分布（可能通过残差）和均值（可能通过合适的初始模型）之间的关系来评估考虑到响应是每株植物的平均产量，假设恒定方差的合理性，当均值较大时，关于均值的变化可能会更大，这似乎是非常合理的。– 格伦_b

每株植物的平均产量严格为正，这意味着我们可以处理它的对数（作为实数）。模型的确定性版本可以有效地重写为：

\ln W_{A} = \ln W_{m A} - b_{A} \ln (1 + a_{A} (N_{A} + α N_{B})) .

$\ln W_{A} = \ln W_{mA} - b_A \ln (1 + a_{A} (N_{A}+\alpha N_{B})).$

明显的随机类比是非线性回归：

\ln W_{A, i} = \ln W_{m A} - b_{A} \ln (1 + a_{A} (N_{A, i} + α N_{B, i})) + ε_{i} ε_{i} \sim IID N (0, σ^{2}) .

$\ln W_{A,i} = \ln W_{mA} - b_A \ln (1 + a_{A} (N_{A,i}+\alpha N_{B,i})) + \varepsilon_i \quad \quad\ \quad \varepsilon_i \sim \text{IID N}(0, \sigma^2).$

这是一个非线性回归，具有未知的系数参数、和，以及未知的误差方差。可以使用以下语法对其进行编程： $a_A$ $b_A$ $\alpha$ $\sigma^2$ R

#Define the formula for the non-linear regression
FUNC    <- function(Wma, Na, Nb, a, b, alpha) { 
                    log(Wma) - b*log(1 + a*(Na + alpha*Nb)) }    
FORMULA <- as.formula(log(Wa) ~ FUNC(Wma, Na, Nb, a, b, alpha));

#Set parameters
a     <- 1.2;
b     <- 0.2;
alpha <- 0.4;
sigma <- 0.1;

#Create mock data for analysis
set.seed(10000);
N   <- 1000;
Wma <- rgamma(N, 20, 4);
Na  <- rgamma(N, 4, 2);
Nb  <- rgamma(N, 6, 1);
Wa  <- rep(0, N);
for (i in 1:N) { Wa[i] <- exp(FUNC(Wma[i], Na[i], Nb[i], a, b, alpha) 
                          + rnorm(1, 0, sigma)) }

DATA <- as.data.frame(cbind(Wa, Wma, Na, Nb));

#Fit the data to the non-linear regression model
MODEL   <- nls(FORMULA, data = DATA, start = list(a = 1, b = 1, alpha = 1));

使用这个特定的模拟数据，模型返回的估计值与用于生成模拟数据的系数值相当接近：

summary(MODEL)
Formula: log(Wa) ~ FUNC(Wma, Na, Nb, a, b, alpha)

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
a      1.30205    0.31614   4.119 4.13e-05 ***
b      0.17763    0.01706  10.410  < 2e-16 ***
alpha  0.53021    0.08063   6.576 7.81e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1028 on 997 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 2.445e-06

其它你可能感兴趣的问题

上一篇展示随时间增长的差异下一篇使用 boostrap 进行偏差和方差估计