基于二维非线性最小二乘(nls)的流行病扩散模型

Question

我试图在R中实现一个流行病扩散模型，扩散的公式是delta_y=(a+b * y) *(N)。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是一个函数，它返回t中的新用户，给定t中的累积用户(实际上是逆累积和函数)。

str(y)传递"int 1:20103，1:16 0 0 0.“。str(delta_y)还提供"int 1:20103，1:16 0 0 0.“。注意，并不是所有的条目都是0，只是最初的许多条目而已。

数据的每个列都有20103个条目。上面的模型适用于单个数据行。

Answer 1

在搜索了Rhelp档案中的错误，并找到了一个邓肯默多克也解决了类似的情况，使用as.vector()将矩阵转换为“长”的as.vector()。，并回顾了在皮涅罗和贝茨的nls和nlsList上的资料，我张贴了一些实验的结果，这些结果可能与您的数据情况相一致。如果我正确理解，您有16种不同的delta_y和y观测“运行”，您希望用同样的非线性模型对它们进行建模。目前还不清楚的是，你是否期望它们都是：(A)具有相同的参数，还是(B)期望系数仅以相同的形式变化。让我们先来看看(A)案件，这就是邓肯默多克( Duncan )9年前提供的“华尔街日报”( as.vector()-solution )所提供的信息。

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

现在在该数据集的所有组上完成，而不考虑group-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