2014-01-14 16 views
0

我無法從存儲在列表中的nls模型的輸出中提取對象。如何從存儲在列表中的nls模型的輸出中提取對象?

這裏是問題:我運行一個漸近模型來估計在(i)年(i)環境中預期要觀察的物種的最大數量。我設法從輸出中提取漸近值而不是p值。

我的數據的一個例子:

env_id year n_sp n_rec 
1 2000 20 20 
1 2000 113 127 
1 2000 170 225 
1 2000 1 1 
1 2000 47 52 
1 2000 6 8 
1 2000 1 1 
1 2000 100 119 
1 2000 30 40 
1 2000 56 60 
1 2000 78 80 
1 2000 34 78 
1 2000 2 2 
1 2000 56 60 
1 2000 89 93 
1 2000 54 67 
1 2000 32 45 
1 2001 7 7 
1 2001 145 162 
1 2001 15 16 
1 2001 24 25 
1 2002 24 25 
1 2002 23 27 
1 2002 128 140 
1 2002 14 14 
1 2002 1 1 
1 2002 1 1 
1 2002 177 189 
1 2002 11 11 
1 2002 3 4 
1 2002 1 1 
1 2002 32 32 
1 2002 10 11 
1 2003 572 834 
1 2003 7 9 
1 2003 4 4 
1 2003 293 396 
1 2003 218 280 
1 2003 12 12 
1 2003 32 33 
1 2003 34 35 
1 2003 10 10 
1 2003 7 7 
1 2003 18 21 
1 2003 2 2 
1 2003 3 3 
1 2003 2 2 
1 2003 74 77 
1 2003 1 1 
2 2000 3 4 
2 2000 6 6 
2 2000 333 470 
2 2000 281 351 
2 2000 4 4 
2 2000 255 319 
2 2000 92 104 
2 2000 218 280 
2 2000 12 12 
2 2000 32 33 
2 2000 34 35 
2 2000 10 10 
2 2000 7 7 
2 2000 18 21 
2 2000 2 2 
2 2000 3 3 
2 2000 2 2 
2 2000 74 77 
2 2000 1 1 
2 2001 88 92 
2 2001 42 44 
2 2001 47 50 
2 2001 5 5 
2 2001 1 1 
2 2001 1 1 
2 2001 1 1 
2 2001 2 2 
2 2001 2 2 
2 2001 2 2 
2 2001 6 7 
2 2001 4 4 
2 2001 13 15 
2 2001 15 16 
2 2003 9 9 
2 2003 94 99 
2 2003 10 10 
2 2003 13 13 
2 2003 2 2 
2 2004 10 10 
2 2004 37 77 
2 2004 23 36 
2 2004 1 1 

運行模型中的代碼如下:

env = unique(dat$env_id) 
res4 = list() 
for (i in 1:length(env)) { 
    dat2 = dat[dat$env_id == env[i],] 
    years2loop <- unique(dat2$year) 
    subres4= list() 
    for (j in 1:length(years2loop)){ 
    subdat2 <- dat2[dat2$year == years2loop [j],] 
    subres4[[j]] = try(nls(n_sp ~ SSasymp(n_rec, Asym, R0, lrc), data = subdat2), silent = T) 
    } 
    for (j in 1:length(subres4)) { 
    if(class(subres4[[j]])=="try-error")subres4[[j]]<-NA 
    } 
    res4[[i]] <- subres4 
} 

這裏是輸出的一個示例:

[[1]] 
[[1]][[1]] 
Nonlinear regression model 
    model: n_sp ~ SSasymp(n_rec, Asym, R0, lrc) 
    data: subdat2 
    Asym  R0  lrc 
577.274 -1.190 -6.445 
residual sum-of-squares: 1476 

Number of iterations to convergence: 3 
Achieved convergence tolerance: 9.236e-07 

[[1]][[2]] 
Nonlinear regression model 
    model: n_sp ~ SSasymp(n_rec, Asym, R0, lrc) 
    data: subdat2 
    Asym  R0  lrc 
1146.9281 0.1631 -7.0899 
residual sum-of-squares: 0.1888 

Number of iterations to convergence: 0 
Achieved convergence tolerance: 6.282e-07 

[[1]][[3]] 
[1] NA 

[[1]][[4]] 
Nonlinear regression model 
    model: n_sp ~ SSasymp(n_rec, Asym, R0, lrc) 
    data: subdat2 
    Asym  R0  lrc 
1841.380 2.166 -7.721 
residual sum-of-squares: 179 

Number of iterations to convergence: 2 
Achieved convergence tolerance: 2.215e-06 


[[2]] 
[[2]][[1]] 
Nonlinear regression model 
    model: n_sp ~ SSasymp(n_rec, Asym, R0, lrc) 
    data: subdat2 
    Asym  R0  lrc 
689.3297 -0.1542 -6.5524 
residual sum-of-squares: 214.6 

Number of iterations to convergence: 1 
Achieved convergence tolerance: 8.802e-06 

[[2]][[2]] 
[1] NA 

[[2]][[3]] 
Nonlinear regression model 
    model: n_sp ~ SSasymp(n_rec, Asym, R0, lrc) 
    data: subdat2 
    Asym  R0  lrc 
819.67028 -0.01942 -6.70025 
residual sum-of-squares: 0.0002441 

Number of iterations to convergence: 0 
Achieved convergence tolerance: 1.348e-06 

[[2]][[4]] 
Nonlinear regression model 
    model: n_sp ~ SSasymp(n_rec, Asym, R0, lrc) 
    data: subdat2 
    Asym  R0  lrc 
48.7765 0.8679 -4.0172 
residual sum-of-squares: 2.614 

Number of iterations to convergence: 0 
Achieved convergence tolerance: 6.642e-06 

我運行下面的代碼來提取漸近值。

rho1<-NULL 
j = 0 
for (i in 1:length(res4)){ 
    for (ii in 1:length(res4[[i]])) { 
    a<-res4[[i]][[ii]] 
    j = j+1 
    if (!is.na(a)) { 
     b<-as.numeric(a$m$getAllPars()[1]) 
     rho1[j]=b 
    } else rho1[j] = NA 
    }} 

這裏我得到的值,但有一個警告消息:

Warning messages: 
1: In if (!is.na(a)) { : 
    the condition has length > 1 and only the first element will be used 

現在,我的問題是如何提取的P值。我想下面的代碼:

p.rho1<-NULL 
for (i in 1:length(res4)){ 
    for (ii in 1:length(res4[[i]])) { 
    c<-res4[[i]][[ii]] 
    j = j+1 
    if (!is.na(c)) { 
     d<-as.numeric(summary(res4[[i]][[ii]])$parameters[1,4]) 
     p.rho1[j]=d 
    } else p.rho1[j] = NA 
    }} 
p.rho1 

有了這個代碼,我得到的警告消息:

Warning messages: 
1: In if (!is.na(c)) { : 
    the condition has length > 1 and only the first element will be used 

此外,我每次我運行代碼時獲得不同數量的值。

有人能幫我解決這個問題嗎?

非常感謝朱莉安娜

+0

在代碼中得到一個單一的列表時, 'sp.accu.t'與'dat'相同?如果我假設並用你的數據運行你的代碼,那麼這些配合都不會收斂。所以我沒有得到適合res4的nls。 – jlhoward

回答

1

您可每次得到一個不同的數值,因爲你沒有在第二組的代碼重新初始化j爲0。

至於警告消息,您試圖辨別該條目是否爲NAnls對象。 nls對象本身是列表,並且is.na應用於列表返回邏輯值的向量,這對應於測試列表中的每個元素。而不是測試不NA,測試是NLS

if(class(c)=="nls") { 

最後,做一些與列表中的每個元素是lapply功能是什麼好

rho1 <- unlist(lapply(res4, function(sublist) { 
    lapply(sublist, function(a) { 
     if(class(a)=="nls") { 
      as.numeric(a$m$getAllPars()[1]) 
     } else { 
      NA 
     } 
    }) 
})) 

p.rho1 <- unlist(lapply(res4, function(sublist) { 
    lapply(sublist, function(a) { 
     if(class(a)=="nls") { 
      as.numeric(summary(a)$parameters[1,4]) 
     } else { 
      NA 
     } 
    }) 
})) 

您也可以使用plyr庫爲更復雜的結構做這種事情。就這樣,的res4創建成爲

library("plyr") 
res4 <- dlply(dat, .(env_id), function(subdat1) { 
    dlply(subdat1, .(year), function(subdat2) { 
     try_default(nls(n_sp ~ SSasymp(n_rec, Asym, R0, lrc), data = subdat2), 
         default = NA, quiet = TRUE) 
    }) 
}) 

如果你並不真的需要列表中,列出了結構,你甚至可以更容易

res4 <- dlply(dat, .(env_id, year), function(subdat2) { 
    try_default(nls(n_sp ~ SSasymp(n_rec, Asym, R0, lrc), data = subdat2), 
       default = NA, quiet = TRUE) 
}) 
相關問題