R：使用二維非線性最小二乘法（nls）的流行擴散模型

Question

我試圖在R中實現流行擴散模型。擴散公式爲delta_y =（a + b * y）*（ NY）。 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是返回噸新用戶給出的累積用戶噸（實際上逆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條目。上述模型適用於單行數據。

Answer 1

在搜索Rhelp存檔後發現該錯誤並找到similar situation was solved by Duncan Murdoch by converting the matrices to "long"-form using as.vector()並在Pinheiro和Bates上查看nls和nlsList上的資料後，我發佈了一些可能與您的數據情況一致的實驗結果。如果我理解正確，你有16個不同的「運行」觀察delta_y和y，你的希望是用相同的非線性模型對它們建模。目前尚不清楚的是你是否期望它們全部：（A）具有相同的參數，或者（B）期望係數僅以相同的形式變化。先來看（A）案例，這是鄧肯默多克9年前提供的解決方案。

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的例子）：

​​

現在在所有組數據集上完成而不考慮組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