R：使用二维非线性最小二乘法（nls）的流行病扩散模型
R: epidemic diffusion model using two-dimensional non-linear least squares (nls)
我正在尝试在 R 中实现流行病扩散模型。扩散的公式是 delta_y=(a+b * y) * (N-y)。 y描述了t时期的当前用户，N描述了潜在用户的数量，delta_y描述了t时期的新用户，a和b是待估计的参数。注意y是前面所有delta_y的累加和。对于单个观察（以 delta_y 和 y 作为向量），模型仅适用于：
model1 <- nls(delta_y ~ (a+b * y) * (N-y))
现在的问题是我有一组这种类型的观察结果，我想为所有这些观察结果估计相同的参数 a 和 b。我试图使用上面的相同公式，但现在 delta_y 和 y 是二维数组而不是向量。我在“qr.qty(QR, resid) 中收到错误：
'qr' 和 'y' 必须具有相同的行数
数据详情：y和delta_y都是16列20103行的二维数组。数组创建如下：
y=matrix(c(data$nearby_1998,data$nearby_1999, data$nearby_2000, ..., data$nearby_2013),nrow=20103)
invCum <- function (data) {result=matrix(nrow=nrow(data), ncol=ncol(data)); result[,1]=data[,1]; for (i in 2:ncol(data)) {result[,i] <- data[,i]-data[,i-1]}; return(result)}
delta_y <- invCum(y)
invCum 是一个函数，returns t 中的新用户给定 t 中的累积用户（实际上是一个逆 cumsum 函数）。
str(y) 交付 "int [1:20103, 1:16] 0 0 0 0 0 0 0 0 0 0 ..."。
str(delta_y) 也提供 "int [1:20103, 1:16] 0 0 0 0 0 0 0 0 0 0 ..."。
请注意，并非所有条目都是 0，只有许多第一个条目。
数据列各有 20103 个条目。上面的模型适用于单行数据。
在 Rhelp 档案中搜索该错误并发现 similar situation was solved by Duncan Murdoch by converting the matrices to "long"-form using as.vector() 并查看 Pinheiro 和 Bates 中 nls 和 nlsList 上的 material 后，我发布了一些实验的结果，这些实验可能与您的数据情况一致。如果我理解正确的话，你有 16 个不同的 "runs" 观察 delta_y 和 y，你希望用相同的非线性模型对它们进行建模。目前尚不清楚的是，您是否期望它们全部：(A) 具有相同的参数，或者 (B) 期望系数仅以相同的形式变化。让我们先来看 (A) 情况，这就是九年前邓肯默多克提供的as.vector()-解决方案所提供的。
newdf <- data.frame( d_y <- as.vector(delta_y),
y = as.vector(y),
grp=rep(letters[1:16], each=20103) )
N= _____ # you need to add this; not sure if it's a constant or vector
# if it varies across groups need to use the rep()-strategy to add to newdf
model1 <- nls( d_y ~ (a+b * y) * (N-y) , data=newdf, start=list(a=0, b=1))
如果另一方面你想要单独的系数：
library(nlme)
model1 <- nlsList(delta_y ~ (a+b * y) * (N-y) | grp, data=newdf, start=c(a=0, b=1) )
这里是一些测试：首先在单个组上（?nls 中的示例）：
DNase1 <- subset(DNase, Run == 1)
> fm2DNase1 <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)),
+ data = DNase1,
+ start = list(xmid = 0, scal = 1) )
> summary(fm2DNase1)
==================
Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal))
Parameters:
Estimate Std. Error t value Pr(>|t|)
xmid -0.02883 0.30785 -0.094 0.927
scal 0.45640 0.27143 1.681 0.115
Residual standard error: 0.3158 on 14 degrees of freedom
Number of iterations to convergence: 14
Achieved convergence tolerance: 1.631e-06
现在在不考虑组 ID 的情况下对该数据集的所有组完成：
> fm2DNase <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)),
+ data = DNase,
+ start = list(xmid = 0, scal = 1) )
> summary(fm2DNase)
==========
Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal))
Parameters:
Estimate Std. Error t value Pr(>|t|)
xmid -0.14816 0.09780 -1.515 0.132
scal 0.46736 0.08691 5.377 2.41e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3291 on 174 degrees of freedom
Number of iterations to convergence: 13
Achieved convergence tolerance: 7.341e-06
最后在每个组上分开，方程形式保持唯一不变：
> fm2DNase <- nlsList(density ~ 1/(1 + exp((xmid - log(conc))/scal))|Run,
+ data = DNase,
+ start = list(xmid = 0, scal = 1) )
> summary(fm2DNase)
Call:
Model: density ~ 1/(1 + exp((xmid - log(conc))/scal)) | Run
Data: DNase
Coefficients:
xmid
Estimate Std. Error t value Pr(>|t|)
10 -0.23467586 0.3527077 -0.66535499 0.4749505
11 -0.18717815 0.3522418 -0.53139112 0.5746396
9 -0.14742434 0.3459987 -0.42608348 0.6521089
1 -0.02882911 0.3403312 -0.08470898 0.9267180
4 -0.01243939 0.3351487 -0.03711604 0.9691708
8 -0.09549007 0.3408348 -0.28016525 0.7741478
5 -0.09216741 0.3367420 -0.27370331 0.7800695
7 -0.25657193 0.3613815 -0.70997535 0.4750054
6 -0.25052019 0.3564816 -0.70275765 0.5051072
2 -0.11218699 0.3245483 -0.34567120 0.7763199
3 -0.23007674 0.3433663 -0.67006203 0.5933597
scal
Estimate Std. Error t value Pr(>|t|)
10 0.4904888 0.3148254 1.557971 0.1076081
11 0.4892928 0.3138277 1.559113 0.1139307
9 0.4723505 0.3075025 1.536087 0.1189793
1 0.4564003 0.3000630 1.521015 0.1148339
4 0.4423467 0.2946883 1.501066 0.1338825
8 0.4582587 0.3018498 1.518168 0.1352101
5 0.4473772 0.2980249 1.501140 0.1407799
7 0.5142468 0.3234251 1.590003 0.1224310
6 0.5007426 0.3185856 1.571768 0.1483103
2 0.4161636 0.2878193 1.445920 0.2457047
3 0.4654567 0.3062277 1.519969 0.2355130
Residual standard error: 0.3491304 on 154 degrees of freedom
我正在尝试在 R 中实现流行病扩散模型。扩散的公式是 delta_y=(a+b * y) * (N-y)。 y描述了t时期的当前用户，N描述了潜在用户的数量，delta_y描述了t时期的新用户，a和b是待估计的参数。注意y是前面所有delta_y的累加和。对于单个观察（以 delta_y 和 y 作为向量），模型仅适用于：
model1 <- nls(delta_y ~ (a+b * y) * (N-y))
现在的问题是我有一组这种类型的观察结果，我想为所有这些观察结果估计相同的参数 a 和 b。我试图使用上面的相同公式，但现在 delta_y 和 y 是二维数组而不是向量。我在“qr.qty(QR, resid) 中收到错误： 'qr' 和 'y' 必须具有相同的行数
数据详情：y和delta_y都是16列20103行的二维数组。数组创建如下：
y=matrix(c(data$nearby_1998,data$nearby_1999, data$nearby_2000, ..., data$nearby_2013),nrow=20103)
invCum <- function (data) {result=matrix(nrow=nrow(data), ncol=ncol(data)); result[,1]=data[,1]; for (i in 2:ncol(data)) {result[,i] <- data[,i]-data[,i-1]}; return(result)}
delta_y <- invCum(y)
invCum 是一个函数，returns t 中的新用户给定 t 中的累积用户（实际上是一个逆 cumsum 函数）。
str(y) 交付 "int [1:20103, 1:16] 0 0 0 0 0 0 0 0 0 0 ..."。 str(delta_y) 也提供 "int [1:20103, 1:16] 0 0 0 0 0 0 0 0 0 0 ..."。 请注意，并非所有条目都是 0，只有许多第一个条目。
数据列各有 20103 个条目。上面的模型适用于单行数据。
在 Rhelp 档案中搜索该错误并发现 similar situation was solved by Duncan Murdoch by converting the matrices to "long"-form using as.vector() 并查看 Pinheiro 和 Bates 中 nls 和 nlsList 上的 material 后，我发布了一些实验的结果，这些实验可能与您的数据情况一致。如果我理解正确的话，你有 16 个不同的 "runs" 观察 delta_y 和 y，你希望用相同的非线性模型对它们进行建模。目前尚不清楚的是，您是否期望它们全部：(A) 具有相同的参数，或者 (B) 期望系数仅以相同的形式变化。让我们先来看 (A) 情况，这就是九年前邓肯默多克提供的as.vector()-解决方案所提供的。
newdf <- data.frame( d_y <- as.vector(delta_y),
y = as.vector(y),
grp=rep(letters[1:16], each=20103) )
N= _____ # you need to add this; not sure if it's a constant or vector
# if it varies across groups need to use the rep()-strategy to add to newdf
model1 <- nls( d_y ~ (a+b * y) * (N-y) , data=newdf, start=list(a=0, b=1))
如果另一方面你想要单独的系数：
library(nlme)
model1 <- nlsList(delta_y ~ (a+b * y) * (N-y) | grp, data=newdf, start=c(a=0, b=1) )
这里是一些测试：首先在单个组上（?nls 中的示例）：
DNase1 <- subset(DNase, Run == 1)
> fm2DNase1 <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)),
+ data = DNase1,
+ start = list(xmid = 0, scal = 1) )
> summary(fm2DNase1)
==================
Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal))
Parameters:
Estimate Std. Error t value Pr(>|t|)
xmid -0.02883 0.30785 -0.094 0.927
scal 0.45640 0.27143 1.681 0.115
Residual standard error: 0.3158 on 14 degrees of freedom
Number of iterations to convergence: 14
Achieved convergence tolerance: 1.631e-06
现在在不考虑组 ID 的情况下对该数据集的所有组完成：
> fm2DNase <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)),
+ data = DNase,
+ start = list(xmid = 0, scal = 1) )
> summary(fm2DNase)
==========
Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal))
Parameters:
Estimate Std. Error t value Pr(>|t|)
xmid -0.14816 0.09780 -1.515 0.132
scal 0.46736 0.08691 5.377 2.41e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3291 on 174 degrees of freedom
Number of iterations to convergence: 13
Achieved convergence tolerance: 7.341e-06
最后在每个组上分开，方程形式保持唯一不变：
> fm2DNase <- nlsList(density ~ 1/(1 + exp((xmid - log(conc))/scal))|Run,
+ data = DNase,
+ start = list(xmid = 0, scal = 1) )
> summary(fm2DNase)
Call:
Model: density ~ 1/(1 + exp((xmid - log(conc))/scal)) | Run
Data: DNase
Coefficients:
xmid
Estimate Std. Error t value Pr(>|t|)
10 -0.23467586 0.3527077 -0.66535499 0.4749505
11 -0.18717815 0.3522418 -0.53139112 0.5746396
9 -0.14742434 0.3459987 -0.42608348 0.6521089
1 -0.02882911 0.3403312 -0.08470898 0.9267180
4 -0.01243939 0.3351487 -0.03711604 0.9691708
8 -0.09549007 0.3408348 -0.28016525 0.7741478
5 -0.09216741 0.3367420 -0.27370331 0.7800695
7 -0.25657193 0.3613815 -0.70997535 0.4750054
6 -0.25052019 0.3564816 -0.70275765 0.5051072
2 -0.11218699 0.3245483 -0.34567120 0.7763199
3 -0.23007674 0.3433663 -0.67006203 0.5933597
scal
Estimate Std. Error t value Pr(>|t|)
10 0.4904888 0.3148254 1.557971 0.1076081
11 0.4892928 0.3138277 1.559113 0.1139307
9 0.4723505 0.3075025 1.536087 0.1189793
1 0.4564003 0.3000630 1.521015 0.1148339
4 0.4423467 0.2946883 1.501066 0.1338825
8 0.4582587 0.3018498 1.518168 0.1352101
5 0.4473772 0.2980249 1.501140 0.1407799
7 0.5142468 0.3234251 1.590003 0.1224310
6 0.5007426 0.3185856 1.571768 0.1483103
2 0.4161636 0.2878193 1.445920 0.2457047
3 0.4654567 0.3062277 1.519969 0.2355130
Residual standard error: 0.3491304 on 154 degrees of freedom