CODEWITHSUNDEEP
r
user2120294
user2120294

Reputation: 3

How to extract objects from an output of an nls model that is stored in a list of lists?

I’m having trouble to extract objects from an output of an nls model that is stored in a list of lists.

Here is the problem: I’m running an asymptotic model to estimate the maximum number of species expected to be observed in environment (i) for the year (j). I managed to extract from the output the asymptotic values but not the p values.

An example of my data:

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

The code for running the model is as follows:

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
}

Here is an example of the output:

[[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

I'm running the following code to extract the asymptotic values.

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
  }}

Here I get the values, but with an warning message:

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

Now, my question is how to extract the p values. I’m trying the following code:

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

With this code I’m getting the warning message: 

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

Additionally, I’m getting a different number of values each time I run the code.

Could somebody help me to fix the problem?

Many thanks, Juliana

Upvotes: 0

Views: 1872

Answers (1)

Brian Diggs
Brian Diggs

Reputation: 58825

You may be getting a different number of values each time because you do not re-initialize j to 0 in the second set of code.

As to the warning message, you are trying to discern whether that entry is NA or an nls object. nls objects are themselves lists, and is.na applied to a list returns a vector of logical values which correspond to testing each element of the list. Rather than testing for not NA, test for is nls

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

Finally, doing something with each element of a list is what the lapply function is good for

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
        }
    })
}))

You can also use the plyr library to do this sort of thing for more complex structures. With that, the creation of res4 becomes

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)
    })
})

If you don't really need a list-of-lists structure, you can get a single list even easier with

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)
})

Upvotes: 2

Related Questions