R Under development (unstable) (2024-06-05 r86695 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## problem reported by Christian Ritz on R-help, Tue, 20 Mar 2007 14:55:38
> 
> ## Dataset PestSci from package drc [w/ 'CURVE' as factor]
> PestSci <-
+ data.frame(CURVE = gl(5, 21),
+            HERBICIDE = rep(c("bentazon", "diuron"), times = c(63, 42)),
+            DOSE = rep(c(0, 0.62, 1.85, 5.56, 16.67, 50, 150, 0, 0.62,
+            1.85, 5.56, 16.67, 50, 150, 0, 0.15, 0.59, 2.34, 9.38, 37.5,
+            150, 0, 0.01, 0.03, 0.1, 0.3, 1, 3, 0, 0.01, 0.03, 0.1, 0.3, 1, 3),
+            each = 3),
+            SLOPE = c(1.81295, 1.86704, 1.95606, 1.39073, 1.15721, 1.06126,
+     0.99409, 0.83298, 0.8334, 0.72513, 0.69548, 0.65299, 0.49855,
+     0.36873, 0.42617, 0.26666, 0.26896, 0.25989, 0.16074, 0.16404,
+     0.1475, 1.02654, 0.91306, 0.89371, 0.59074, 0.669, 0.5965,
+     0.37561, 0.44823, 0.42093, 0.31874, 0.27342, 0.2725, 0.27182,
+     0.21752, 0.19981, 0.17332, 0.17949, 0.15623, 0.12855, 0.14524,
+     0.11533, 1.03872, 1.0917, 1.10324, 0.94274, 0.91256, 1.02352,
+     0.78689, 0.69706, 0.65989, 0.5372, 0.51324, 0.54981, 0.37401,
+     0.34033, 0.32491, 0.30518, 0.24593, 0.289, 0.17414, 0.12275,
+     0.14788, 2.20963, 2.27931, 2.14703, 2.18831, 2.08863, 2.06676,
+     2.18827, 2.10748, 1.84474, 1.78805, 1.75547, 1.61381, 0.70295,
+     0.6983, 0.74045, 0.20673, 0.20784, 0.22402, 0.05268, 0.06519,
+     0.09258, 1.94033, 1.80193, 1.71586, 1.71586, 1.98471, 1.74905,
+     1.87795, 1.64081, 1.53094, 1.50709, 1.41035, 1.35367, 0.64427,
+     0.62185, 0.60337, 0.14073, 0.12928, 0.15016, 8e-05, 0.00262,
+     0.00303))
> 
> 
> library(nlme)
> sv <- c(0.43355869, 2.49963220, 0.05861799, 1.73290589, 0.38153146, 0.24316978)
> 
> ## bug only triggered from the  " factor(.) " below:
> fm <- nlme(SLOPE ~ c + (d-c)/(1+exp(b*(log(DOSE)-log(e)))),
+            fixed = list(b ~ factor(HERBICIDE)-1,
+                         c ~ 1,
+                         d ~ 1,
+                         e ~ factor(HERBICIDE)-1),
+            random = d ~ 1 | CURVE,
+            start = sv, data = PestSci)
> fm
Nonlinear mixed-effects model fit by maximum likelihood
  Model: SLOPE ~ c + (d - c)/(1 + exp(b * (log(DOSE) - log(e)))) 
  Data: PestSci 
  Log-likelihood: 113.6823
  Fixed: list(b ~ factor(HERBICIDE) - 1, c ~ 1, d ~ 1, e ~ factor(HERBICIDE) -      1) 
b.factor(HERBICIDE)bentazon   b.factor(HERBICIDE)diuron 
                 0.56353864                  1.79126258 
                          c                           d 
                 0.05198172                  1.57545308 
e.factor(HERBICIDE)bentazon   e.factor(HERBICIDE)diuron 
                 1.57142104                  0.20111247 

Random effects:
 Formula: d ~ 1 | CURVE
                d   Residual
StdDev: 0.4699677 0.07139796

Number of Observations: 105
Number of Groups: 5 
> ## failed in contrMat / contr.treatment in 3.1-78
> 
> stopifnot(
+     all.equal(as.numeric(VarCorr(fm)[,"StdDev"]),
+ 	      c(0.46996772, 0.07139796), tol = 5e-5)
+     ,
+     all.equal(as.numeric(fixef(fm)),
+ 	      c(0.563538644, 1.7912626,
+ 		0.0519817224, 1.57545308,
+ 		1.5714210, 0.201112468), tol = 1e-7)
+ )
> (sfm <- summary(fm))
Nonlinear mixed-effects model fit by maximum likelihood
  Model: SLOPE ~ c + (d - c)/(1 + exp(b * (log(DOSE) - log(e)))) 
  Data: PestSci 
        AIC       BIC   logLik
  -211.3646 -190.1329 113.6823

Random effects:
 Formula: d ~ 1 | CURVE
                d   Residual
StdDev: 0.4699677 0.07139796

Fixed effects:  list(b ~ factor(HERBICIDE) - 1, c ~ 1, d ~ 1, e ~ factor(HERBICIDE) -      1) 
                                Value  Std.Error DF   t-value p-value
b.factor(HERBICIDE)bentazon 0.5635386 0.03631152 95 15.519554  0.0000
b.factor(HERBICIDE)diuron   1.7912626 0.10634500 95 16.843882  0.0000
c                           0.0519817 0.02210700 95  2.351369  0.0208
d                           1.5754531 0.21706624 95  7.257937  0.0000
e.factor(HERBICIDE)bentazon 1.5714210 0.19138896 95  8.210615  0.0000
e.factor(HERBICIDE)diuron   0.2011125 0.00771419 95 26.070469  0.0000
 Correlation: 
                            b.fctr(HERBICIDE)b b.fctr(HERBICIDE)d c      d     
b.factor(HERBICIDE)diuron    0.274                                             
c                            0.601              0.456                          
d                           -0.021             -0.019             -0.011       
e.factor(HERBICIDE)bentazon  0.025             -0.221             -0.486 -0.032
e.factor(HERBICIDE)diuron   -0.256             -0.013             -0.426 -0.011
                            e.fctr(HERBICIDE)b
b.factor(HERBICIDE)diuron                     
c                                             
d                                             
e.factor(HERBICIDE)bentazon                   
e.factor(HERBICIDE)diuron    0.208            

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.40381491 -0.54615415 -0.06518212  0.56616131  2.73251640 

Number of Observations: 105
Number of Groups: 5 
> cfT <- sfm$tTable
> stopifnot(is.matrix(cfT), identical(cfT[,"Value"], fixef(fm)),
+ 	  all.equal(as.vector(cfT[,"Std.Error"]),
+ 		    c(0.0363115, 0.106345, 0.022107,
+ 		      0.217066, 0.191389, 0.00771419), tol = 5e-5)
+ 	  ,
+ 	  all.equal(as.numeric(sfm[c("AIC", "BIC", "logLik")]),
+ 		    c(-211.36462, -190.13294, 113.68231), tol = 1e-7)
+ 	  )
> 
> proc.time()
   user  system elapsed 
   0.20    0.04    0.25