Skip to content

Error generated by nlme::gnls

9 messages · Aaron Crowley, Ben Tupper, Rolf Turner +3 more

Original
Aaron Crowley
# 
Approximately 6 months ago, I successfully performed a model-based
network meta-analysis using the gnls function from the nlme package in
R (nlme::gnls). In this analysis, a binary response is captured as a
non-linear non-parametric placebo response model combined with drug
effect as a function of dose and time and is fitted on aggregate data
from 57 trials and 13 treatments.

I am now trying to rerun the analysis but the code no longer runs. An
error message is generated when trying to fit the model using
nlme::gnls. A reproducible example generating the same error message
is below.

library(nlme)

set.seed(548)

df <- Soybean
df$x01 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x02 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x03 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x04 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x05 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x06 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x07 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x08 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x09 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x10 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x11 <- sample(c(0, 1), size = nrow(df), replace = TRUE)
df$x12 <- sample(c(0, 1), size = nrow(df), replace = TRUE)

gnls(
  weight ~ x,
  data = df,
  params = (x ~ -1 + x01 + x02 + x03 + x04 + x05 + x06 + x07 + x08 +
x09 + x10 + x11 + x12),
  start = rep(0, 13)
)
# Error in if (deparse(params[[nm]][[3]]) != "1") { :
#     the condition has length > 1

sessionInfo()
# R version 4.2.1 (2022-06-23 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 22000)
#
# Matrix products: default
#
# locale:
#   [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United
States.utf8
# [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
# [5] LC_TIME=English_United States.utf8
#
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base
#
# other attached packages:
#   [1] nlme_3.1-157
#
# loaded via a namespace (and not attached):
#   [1] compiler_4.2.1  tools_4.2.1     grid_4.2.1      lattice_0.20-45

Aaron Crowley
Principal Scientist, Biostatistics

e.  aaron at genesisrg.com
w. genesisrg.com
2 days later
Ben Tupper
# 
Could this be related to a new if() behavior introduced in v4.2.0 ?
See the "SIGNIFICANT USER-VISIBLE CHANGES" for v4.2.0 in the NEWS

https://cloud.r-project.org/doc/manuals/r-release/NEWS.html
On Sat, Jul 23, 2022 at 6:26 PM Aaron Crowley <aaron at genesisrg.com> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rolf Turner
    

    

      
      #
    

  


  

      
On Sat, 23 Jul 2022 21:00:25 -0400

        
Ben Tupper <btupper at bigelow.org> wrote:

        

            

      
      
      
    

        
No.  What's going on is much weirder than that, and looks to me like
a bug has been introduced in formula() or in something related.

My impression is that if the right hand side of a formula gets "too
long", then it gets split into parts --- which messes everything up.

Consider:

lhs1 <- paste(paste0("x",1:12),collapse=" + ")
f1   <- as.formula(paste("x ~ -1 +",lhs1))
print(length(deparse(f1[[3]])))

You get the value 2.  Note that deparse(f1[[3]]) gives

[1] "-1 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + "
[2] "    x12"       

If we shorten the right hand side of the formula just a wee bit

lhs2 <- paste(paste0("x",1:11),collapse=" + ")
f2   <- as.formula(paste("x ~ -1 +",lhs2))

then deparse(f2[[3]]) is of length 1, as it should be:

            

      
      
      
    

        
I could of course be missing something, but it really looks to me as if
something has gone up to Puttee here.

Some input from someone in R-Core would be valuable here.

cheers,

Rolf Turner

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ivan Krylov
    

    

      
      #
    

  


  

      
On Sun, 24 Jul 2022 16:03:24 +1200

        
Rolf Turner <r.turner at auckland.ac.nz> wrote:

        

            

      
      
      
    

        
For new enough R (? 4.0), it's possible to use deparse1() [*], which
guarantees to return a single string. Otherwise, the simplest
workaround would be paste(deparse(...), collapse = ' '), which is
similar to how deparse1 is defined.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rolf Turner
    

    

      
      #
    

  


  

      
On Sun, 24 Jul 2022 07:58:01 +0300

        
Ivan Krylov <krylov.r00t at gmail.com> wrote:

        

            

      
      
      
    

        
Yes but .....

I guess my posting did not make things clear.  Sorry 'bout that! My
posting, and Ben Tupper's posting, were in response to an earlier
posting by Aaron Crowley which dealt with an error that now results
from a call to gnls(). Apparently the exact same call had worked
without error in the past.

The call to deparse(), which triggers the error, is in the code
of gnls().  (Line 223 of gnls.R.)

The maintainer of the nlme package (who is, according to maintainer(),
"R-core") could change the code so that it uses invokes deparse1()
rather than deparse, but the user cannot do so, not without in effect
re-creating the package.

Also, the question remains:  why did Aaron Crowley's code work in the
past, whereas now it throws an error?  What changed?

cheers,

Rolf Turner

P.S. Aaron Crowley tells us that he is currently running R 4.2.1.
However he did not tell us what version of R he was running at the time
("approximately six months ago") when the call to gnls() executed
successfully.  Perhaps that version was older than 4.0.  Even so, it's
all a bit mysterious.

R. T.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ivan Krylov
    

    

      
      #
    

  


  

      
Sorry for being too terse in my previous e-mail!

On Sun, 24 Jul 2022 23:03:02 +1200

        
Rolf Turner <r.turner at auckland.ac.nz> wrote:

        

            

      
      
      
    

        
You're right. I think there's a buglet in nlme::gnls that nobody
noticed until R 4.2.0 was released *and* Aaron Crowley used the
function with a sufficiently long formula.

            

      
      
      
    

        
gnls() may have been performing the `if (deparse(...) != '1')` test for
a long time, but never crashed before because it wasn't a fatal error
until R 4.2.0. Previously, if() would issue a warning and use the first
element of the boolean vector.

R 4.2.0 was released this April, which was less than 6 months ago. I
think it all fits.

A temporary solution would be to make use of the fact that R is a very
dynamic language and perform surgery on a live function inside a loaded
package:

library(codetools)

nlme <- loadNamespace('nlme')
unlockBinding('gnls', nlme)
nlme$gnls <- `body<-`(fun = nlme$gnls, value = walkCode(
	body(nlme$gnls), makeCodeWalker(
		call = function(e, w)
			as.call(lapply(as.list(e), function(ee)
				if (!missing(ee)) walkCode(ee, w)
			)),
		leaf = function(e, w)
			if (is.symbol(e) && e == 'deparse') {
				as.name('deparse1')
			} else e
	)
))
lockBinding('gnls', nlme)
rm(nlme)

grep('deparse', deparse(nlme::gnls), value = TRUE)
# [1] "                deparse1(pp[[3]]), sep = \"~\"), collapse =
\",\"), " # [2] "        if (deparse1(params[[nm]][[3]]) != \"1\") {"
# [3] "        list(row.names(dataModShrunk), deparse1(form[[2]]))), "

Aaron's example seems to work after this, modulo needing 12 starting
values instead of 13.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Bill Dunlap
    

    

      
      #
    

  


  

      
I think the intent of this code was to see if the formula had solely a
literal 1 on the right hand side.
Then !identical(pp[[3]], 1) would do it, avoiding the overhead of calling
deparse.  Note that the
1 might be an integer, which the current code would (erroneously) not
catch, so
   !identical(pp[[3]], 1) && !identical(pp[[3]], 1L)
or something equivalent would be better.

-Bill

        
On Sun, Jul 24, 2022 at 5:58 AM Ivan Krylov <krylov.r00t at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Rolf Turner
    

    

      
      #
    

  


  

      
On Sun, 24 Jul 2022 15:57:25 +0300

        
Ivan Krylov <krylov.r00t at gmail.com> wrote:

        

            

      
      
      
    

        
Makes sense. :-)

            

      
      
      
    

        
N'ya-ha!  Makes sense too.

            

      
      
      
    

        
Indeed,

            

      
      
      
    

        
Such surgery is beyond the capability of us ordinary mortals, but given
your explicit and clear recipe it should be do-able. 

Thanks for all of this.

cheers,

Rolf

P. S.  Ben:  you were correct in your original conjecture, to which I
erroneously said "No".  Lack of insight on my part.

R.

  
    

      
      
    

  


  


      

    

    
3 days later

    

      
      

        


  

        

      Martin Maechler
    

    

      
      #
    

  


  

          

      
      
      
    

        
> I think the intent of this code was to see if the formula
    > had solely a literal 1 on the right hand side.  Then
    > !identical(pp[[3]], 1) would do it, avoiding the overhead
    > of calling deparse.  Note that the 1 might be an integer,
    > which the current code would (erroneously) not catch, so
    > !identical(pp[[3]], 1) && !identical(pp[[3]], 1L) or
    > something equivalent would be better.

    > -Bill

Thanks a lot, Bill, Ivan, Ben  et al. !

This is indeed a very old coding bug triggered by the more
strict checks in  R 4.2.x.

I will indeed try Bill's proposal rather than remaining with
deparse by using deparse1().

"Of course", this should hopefully be fixed in the next release
of nlme.

Martin Maechler
ETH Zurich   and  R Core team


    > On Sun, Jul 24, 2022 at 5:58 AM Ivan Krylov

        
> <krylov.r00t at gmail.com> wrote:

        
>> Sorry for being too terse in my previous e-mail!
    >> 
    >> On Sun, 24 Jul 2022 23:03:02 +1200 Rolf Turner

        
>> <r.turner at auckland.ac.nz> wrote:

        
>> 
    >> > The maintainer of the nlme package (who is, according
    >> to maintainer(), > "R-core") could change the code so
    >> that it uses invokes deparse1() > rather than deparse,
    >> but the user cannot do so, not without in effect >
    >> re-creating the package.
    >> 
    >> You're right. I think there's a buglet in nlme::gnls that
    >> nobody noticed until R 4.2.0 was released *and* Aaron
    >> Crowley used the function with a sufficiently long
    >> formula.
    >> 
    >> > Also, the question remains: why did Aaron Crowley's
    >> code work in the > past, whereas now it throws an error?
    >> What changed?
    >> 
    >> gnls() may have been performing the `if (deparse(...) !=
    >> '1')` test for a long time, but never crashed before
    >> because it wasn't a fatal error until R
    >> 4.2.0. Previously, if() would issue a warning and use the
    >> first element of the boolean vector.
    >> 
    >> R 4.2.0 was released this April, which was less than 6
    >> months ago. I think it all fits.
    >> 
    >> A temporary solution would be to make use of the fact
    >> that R is a very dynamic language and perform surgery on
    >> a live function inside a loaded package:
    >> 
    >> library(codetools)
    >> 
    >> nlme <- loadNamespace('nlme') unlockBinding('gnls', nlme)
    >> nlme$gnls <- `body<-`(fun = nlme$gnls, value = walkCode(
    >> body(nlme$gnls), makeCodeWalker( call = function(e, w)
    >> as.call(lapply(as.list(e), function(ee) if (!missing(ee))
    >> walkCode(ee, w) )), leaf = function(e, w) if
    >> (is.symbol(e) && e == 'deparse') { as.name('deparse1') }
    >> else e ) )) lockBinding('gnls', nlme) rm(nlme)
    >> 
    >> grep('deparse', deparse(nlme::gnls), value = TRUE) # [1]
    >> " deparse1(pp[[3]]), sep = \"~\"), collapse = \",\"), " #
    >> [2] " if (deparse1(params[[nm]][[3]]) != \"1\") {" # [3]
    >> " list(row.names(dataModShrunk), deparse1(form[[2]]))), "
    >> 
    >> Aaron's example seems to work after this, modulo needing
    >> 12 starting values instead of 13.
    >> 
    >> --
    >> Best regards, Ivan
    >> 
    >> ______________________________________________
    >> R-help at r-project.org mailing list -- To UNSUBSCRIBE and
    >> more, see https://stat.ethz.ch/mailman/listinfo/r-help
    >> PLEASE do read the posting guide
    >> http://www.R-project.org/posting-guide.html and provide
    >> commented, minimal, self-contained, reproducible code.
    >> 

    > 	[[alternative HTML version deleted]]

    > ______________________________________________
    > R-help at r-project.org mailing list -- To UNSUBSCRIBE and
    > more, see https://stat.ethz.ch/mailman/listinfo/r-help
    > PLEASE do read the posting guide
    > http://www.R-project.org/posting-guide.html and provide
    > commented, minimal, self-contained, reproducible code.