Benjamin Christoffersen
2018-May-05 10:32 UTC
[Rd] Bug in profile.nls with algorithm = "plinear"
Dear sirs It seems like there is a bug in `profile.nls` with `algorithm "plinear"` when a matrix is supplied on the right hand side. Here is the bug and a potential fix ##### # example where profile.nls does not work with `plinear` but does with # `default` require(graphics) set.seed(1) DNase1 <- subset(DNase, Run == 1) x <- rnorm(nrow(DNase1)) f1 <- nls(density ~ b1/(1 + exp((xmid - log(conc))/scal)) + b2 / x, data = DNase1, start = list(xmid = 0, scal = 1, b1 = 1, b2 = 1)) coef (f1) #R> xmid scal b1 b2 #R> 1.5055308 1.0455951 2.3707962 0.0006887 confint(f1) #R> Waiting for profiling to be done... #R> 2.5% 97.5% #R> xmid 1.323034 1.727130 #R> scal 0.974952 1.123036 #R> b1 2.193703 2.600172 #R> b2 -0.001597 0.002978 f2 <- nls(density ~ cbind(1/(1 + exp((xmid - log(conc))/scal)), x), data = DNase1, start = list(xmid = 0, scal = 1), algorithm = "plinear") coef (f2) #R> xmid scal .lin1 .lin.x #R> 1.461636 1.028726 2.323707 0.008807 confint(f2) # this fails #R> Waiting for profiling to be done... #R> Error in attr(ans, "gradient")[c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, : #R> (subscript) logical subscript too long traceback() # [output output abbreviated] ##### # It seems to fail due to the dimension in gradCall environment(f1$m$setPars)$gradCall #R> attr(ans, "gradient")[c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, #R> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE), #R> c(TRUE, TRUE, TRUE, TRUE), drop = FALSE] (ca <- environment(f2$m$setPars)$gradCall) #R> attr(ans, "gradient")[c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, #R> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE), #R> c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, #R> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE), c(TRUE, TRUE)] f2$m$setVarying # [output output abbreviated] environment(f2$m$setVarying)$getRHS.varying #R> function () #R> { #R> ans <- getRHS.noVarying() #R> attr(ans, "gradient") <- eval(gradCall) #R> ans #R> } #R> <bytecode: 0x000000000ce31768> #R> <environment: 0x000000000cbdbad0> dim(attr(environment(f2$m$setVarying)$getRHS.noVarying(), "gradient")) #R> [1] 16 2 2 sapply(ca[3:5], length) #R> [1] 16 16 2 ##### # It seems like the `nlsModel.plinear` in `stats/R/nls.R` should be # # # ... # for(i in 1:marg + 1L) # gradSetArgs[[i]] <- rep_len(TRUE, dimGrad[i-1]) # # ... # gradCall <- # switch(length(gradSetArgs) - 1L, # call("[", gradSetArgs[[1L]], drop = FALSE, gradSetArgs[[2L]]), # call("[", gradSetArgs[[1L]], drop = FALSE, gradSetArgs[[2L]], gradSetArgs[[3L]]), # call("[", gradSetArgs[[1L]], drop = FALSE, gradSetArgs[[2L]], gradSetArgs[[3L]], # gradSetArgs[[4L]]), # call("[", gradSetArgs[[1L]], drop = FALSE, gradSetArgs[[2L]], gradSetArgs[[3L]], # gradSetArgs[[4L]], gradSetArgs[[5L]])) # this seems to fix the problem. source("R/nls_mod.R") # `nls_mod` is assigned in `R/nls_mod.R` f3 <- nls_mod(density ~ cbind(1/(1 + exp((xmid - log(conc))/scal)), x), data = DNase1, start = list(xmid = 0, scal = 1), algorithm = "plinear") confint(f3) #R> Waiting for profiling to be done... #R> 2.5% 97.5% #R> xmid 1.3169 1.635 #R> scal 0.9667 1.096 ##### # The example on the help page still gives the same ca <- quote(nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)), data = DNase1, start = list(xmid = 0, scal = 1), algorithm = "plinear")) confint(eval(ca)) #R> Waiting for profiling to be done... #R> 2.5% 97.5% #R> xmid 1.3225 1.677 #R> scal 0.9748 1.114 ca[[1]] <- quote(nls_mod) confint(eval(ca)) #R> Waiting for profiling to be done... #R> 2.5% 97.5% #R> xmid 1.3225 1.677 #R> scal 0.9748 1.114 ##### # session info sessionInfo() #R> R version 3.5.0 (2018-04-23) #R> Platform: x86_64-w64-mingw32/x64 (64-bit) #R> Running under: Windows >= 8 x64 (build 9200) #R> #R> Matrix products: default #R> #R> locale: #R> [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 #R> [4] LC_NUMERIC=C LC_TIME=English_United States.1252 #R> #R> attached base packages: #R> [1] stats graphics grDevices utils datasets methods base #R> #R> loaded via a namespace (and not attached): #R> [1] MASS_7.3-49 compiler_3.5.0 tools_3.5.0 packrat_0.4.8-1