【问题标题】:How to extract objects from an output of an nls model that is stored in a list of lists?如何从存储在列表列表中的 nls 模型的输出中提取对象?
【发布时间】:2014-11-29 05:53:31
【问题描述】:

我无法从存储在列表列表中的 nls 模型的输出中提取对象。

问题出在这里:我正在运行一个渐近模型来估计预计在环境 (i) 中 (j) 年可观察到的最大物种数量。我设法从输出中提取渐近值而不是 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

此外,每次运行代码时,我都会得到不同数量的值。

有人可以帮我解决这个问题吗?

非常感谢,朱莉安娜

【问题讨论】:

  • 在您的代码中,sp.accu.t 是否与dat 相同?如果我假设并使用您的数据运行您的代码,那么所有拟合都不会收敛。所以我没有得到适合 res4 的 nls。

标签: r


【解决方案1】:

由于您没有在第二组代码中将 j 重新初始化为 0,因此您每次可能会获得不同数量的值。

关于警告消息,您正试图辨别该条目是NA 还是nls 对象。 nls 对象本身就是列表,is.na 应用于列表会返回一个逻辑值向量,该向量对应于测试列表中的每个元素。与其测试 not NA,不如测试 is 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)
})

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2016-07-12
    • 1970-01-01
    • 2020-03-17
    • 1970-01-01
    相关资源
    最近更新 更多