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. > #### Testing Examples --- related to the new "sigma = fixed" feature > #### ================================================================ > > ## library(nlme,lib.loc=lib.loc) > library(nlme) > > nlme:::doExtras() ## used below {when activated by the tester, e.g., MM..r} [1] FALSE > ## isSun <- Sys.info()[["sysname"]] == "SunOS" > > ##=== example 1 general linear model page 251 gls ML and LME ================ > ## > ex <- "ex1_gls_page251"; .pt <- proc.time() > ## > cat("\n example ", ex,"\n") example ex1_gls_page251 > sigma <- 2 > cat("\nFixed sigma= ",sigma," estimation method 'ML'\n") Fixed sigma= 2 estimation method 'ML' > t1.fix.ML.gls <- gls(distance ~ Sex *I(age-11), data = Orthodont, + correlation = corSymm(form = ~1 | Subject), + weights = varIdent(form = ~1 |age), + control = glsControl(sigma = sigma), + method = "ML") > (s1M <- summary(t1.fix.ML.gls)) Generalized least squares fit by maximum likelihood Model: distance ~ Sex * I(age - 11) Data: Orthodont AIC BIC logLik 446.342 481.2097 -210.171 Correlation Structure: General Formula: ~1 | Subject Parameter estimate(s): Correlation: 1 2 3 2 0.499 3 0.607 0.526 4 0.473 0.701 0.708 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | age Parameter estimates: 8 10 12 14 1.0000000 0.9624373 1.1663059 1.0431039 Coefficients: Value Std.Error t-value p-value (Intercept) 25.000370 0.4326141 57.78908 0.0000 SexFemale -2.342983 0.6777760 -3.45687 0.0008 I(age - 11) 0.822023 0.0796287 10.32320 0.0000 SexFemale:I(age - 11) -0.345054 0.1247543 -2.76587 0.0067 Correlation: (Intr) SexFml I(-11) SexFemale -0.638 I(age - 11) 0.183 -0.117 SexFemale:I(age - 11) -0.117 0.183 -0.638 Standardized residuals: Min Q1 Med Q3 Max -2.76715021 -0.70292159 -0.08541963 0.67916836 2.48284979 Residual standard error: 2 Degrees of freedom: 108 total; 104 residual > (a1M <- anova(t1.fix.ML.gls)) # sequential Denom. DF: 104 numDF F-value p-value (Intercept) 1 5094.398 <.0001 Sex 1 9.007 0.0034 I(age - 11) 1 123.585 <.0001 Sex:I(age - 11) 1 7.650 0.0067 > (a1Mm<- anova(t1.fix.ML.gls, type = "marginal")) Denom. DF: 104 numDF F-value p-value (Intercept) 1 3339.578 <.0001 Sex 1 11.950 0.0008 I(age - 11) 1 106.568 <.0001 Sex:I(age - 11) 1 7.650 0.0067 > ## t_{n} ^2 == F_{1,n}: > stopifnot(all.equal(as.vector(s1M$tTable[,"t-value"] ^ 2), + a1Mm[,"F-value"], tolerance = 1e-14), + identical(2, sigma(t1.fix.ML.gls))) > > ## > cat("\nFixed sigma= ",sigma," estimation method 'REML'\n") Fixed sigma= 2 estimation method 'REML' > t1.fix.REML.gls <- gls(distance ~ Sex*I(age-11), data = Orthodont, + correlation = corSymm(form = ~1 | Subject), + weights = varIdent(form = ~1 |age), + control = glsControl(sigma = sigma, + maxIter = 1000, msMaxIter = 200), + method = "REML") > (s1R <- summary(t1.fix.REML.gls)) Generalized least squares fit by REML Model: distance ~ Sex * I(age - 11) Data: Orthodont AIC BIC logLik 451.8351 486.2122 -212.9176 Correlation Structure: General Formula: ~1 | Subject Parameter estimate(s): Correlation: 1 2 3 2 0.509 3 0.603 0.536 4 0.465 0.701 0.716 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | age Parameter estimates: 8 10 12 14 1.0000000 0.9804346 1.1796126 1.0757388 Coefficients: Value Std.Error t-value p-value (Intercept) 25.012045 0.4324472 57.83838 0.0000 SexFemale -2.356145 0.6775145 -3.47763 0.0007 I(age - 11) 0.821141 0.0809370 10.14543 0.0000 SexFemale:I(age - 11) -0.344059 0.1268039 -2.71332 0.0078 Correlation: (Intr) SexFml I(-11) SexFemale -0.638 I(age - 11) 0.204 -0.130 SexFemale:I(age - 11) -0.130 0.204 -0.638 Standardized residuals: Min Q1 Med Q3 Max -2.77431180 -0.68579222 -0.08434759 0.65942065 2.47568820 Residual standard error: 2 Degrees of freedom: 108 total; 104 residual > (a1R <- anova(t1.fix.REML.gls)) Denom. DF: 104 numDF F-value p-value (Intercept) 1 5115.788 <.0001 Sex 1 8.923 0.0035 I(age - 11) 1 119.455 <.0001 Sex:I(age - 11) 1 7.362 0.0078 > (a1Rm <- anova(t1.fix.REML.gls, type="marginal")) Denom. DF: 104 numDF F-value p-value (Intercept) 1 3345.278 <.0001 Sex 1 12.094 0.0007 I(age - 11) 1 102.930 <.0001 Sex:I(age - 11) 1 7.362 0.0078 > intervals(t1.fix.REML.gls) # now work (via 'apVar') Approximate 95% confidence intervals Coefficients: lower est. upper (Intercept) 24.1544861 25.0120451 25.86960414 SexFemale -3.6996814 -2.3561446 -1.01260786 I(age - 11) 0.6606394 0.8211405 0.98164158 SexFemale:I(age - 11) -0.5955162 -0.3440593 -0.09260248 Correlation structure: lower est. upper cor(1,2) 0.2987693 0.5086627 0.6716156 cor(1,3) 0.4126106 0.6025511 0.7422774 cor(1,4) 0.2482924 0.4648599 0.6371695 cor(2,3) 0.3210549 0.5364025 0.6990086 cor(2,4) 0.5141701 0.7014643 0.8249036 cor(3,4) 0.5385465 0.7161527 0.8328285 Variance function: lower est. upper 10 0.7685047 0.9804346 1.250808 12 0.9421573 1.1796126 1.476915 14 0.8454639 1.0757388 1.368733 Residual standard error: [1] 2 2 2 > ## t_{n} ^2 == F_{1,n}: > stopifnot(all.equal(as.vector(s1R$tTable[,"t-value"] ^ 2), + a1Rm[,"F-value"], tolerance = 1e-14)) > > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.39 0.01 0.41 > > ##=== example 2 linear mixed model page 147 lme ML and REML ================ > ex <- "ex2_lme_page147"; .pt <- proc.time() > ## > cat("\n example ", ex,"\n") example ex2_lme_page147 > ## > method <- "ML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method ML > t1.fix.ML.lme <- lme(distance ~ I(age-11), data = Orthodont, + control = lmeControl(sigma = sigma), + method = method) > summary (t1.fix.ML.lme) Linear mixed-effects model fit by maximum likelihood Data: Orthodont AIC BIC logLik 458.7204 472.1311 -224.3602 Random effects: Formula: ~I(age - 11) | Subject Structure: General positive-definite StdDev Corr (Intercept) 2.1330288 (Intr) I(age - 11) 0.2863615 0.383 Residual 1.0000000 Fixed effects: distance ~ I(age - 11) Value Std.Error DF t-value p-value (Intercept) 24.023148 0.4255878 80 56.44699 0 I(age - 11) 0.660185 0.0705779 80 9.35399 0 Correlation: (Intr) I(age - 11) 0.294 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.72724194 -0.52002269 0.04191596 0.50737507 5.03926723 Number of Observations: 108 Number of Groups: 27 > anova (t1.fix.ML.lme) numDF denDF F-value p-value (Intercept) 1 80 3156.0348 <.0001 I(age - 11) 1 80 87.4972 <.0001 > intervals(t1.fix.ML.lme) Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 23.1840802 24.0231481 24.8622161 I(age - 11) 0.5210372 0.6601852 0.7993331 Random Effects: Level: Subject lower est. upper sd((Intercept)) 1.67885287 2.1330288 2.7100717 sd(I(age - 11)) 0.19442377 0.2863615 0.4217740 cor((Intercept),I(age - 11)) -0.02108675 0.3829427 0.6794528 > stopifnot(sigma(t1.fix.ML.lme) == 1) > > method <- "REML" > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method REML > t1.fix.REML.lme <- lme(distance ~ I(age-11), data = Orthodont, + control = lmeControl(sigma = sigma), method = method) > summary (t1.fix.REML.lme) Linear mixed-effects model fit by REML Data: Orthodont AIC BIC logLik 462.2497 475.5669 -226.1249 Random effects: Formula: ~I(age - 11) | Subject Structure: General positive-definite StdDev Corr (Intercept) 2.1600597 (Intr) I(age - 11) 0.2918805 0.38 Residual 1.0000000 Fixed effects: distance ~ I(age - 11) Value Std.Error DF t-value p-value (Intercept) 24.023148 0.4266952 80 56.30049 0 I(age - 11) 0.660185 0.0707615 80 9.32972 0 Correlation: (Intr) I(age - 11) 0.294 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.72541865 -0.52604171 0.04033649 0.50851493 5.03597974 Number of Observations: 108 Number of Groups: 27 > anova (t1.fix.REML.lme) numDF denDF F-value p-value (Intercept) 1 80 3139.6736 <.0001 I(age - 11) 1 80 87.0437 <.0001 > intervals(t1.fix.REML.lme) Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 23.1739976 24.0231481 24.8722987 I(age - 11) 0.5193653 0.6601852 0.8010051 Random Effects: Level: Subject lower est. upper sd((Intercept)) 1.69468705 2.1600597 2.7532269 sd(I(age - 11)) 0.19819342 0.2918805 0.4298539 cor((Intercept),I(age - 11)) -0.02517111 0.3799697 0.6779037 > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.07 0 0.08 > > ##=== example 3 general non-linear model page 402/ page 512 gnls ls ======== > ex <- "ex3_gnls_page402"; .pt <- proc.time() > ## > cat("\n example ", ex,"\n") example ex3_gnls_page402 > ## > method <- "LS" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method LS > t1.fix.gnls <- gnls( rate ~SSasympOff(pressure, Asym, lrc, c0), data = Dialyzer, + params = list(Asym + lrc ~ QB, c0 ~ 1), + start = c(53.6,8.6,0.51,-0.26, 0.225), + control = gnlsControl(sigma = 1)) > stopifnot(is.null(t1.fix.gnls$apVar)) ## as is has *no* varying ranef-parameters > summary (t1.fix.gnls) Generalized nonlinear least squares fit Model: rate ~ SSasympOff(pressure, Asym, lrc, c0) Data: Dialyzer AIC BIC logLik 2206.684 2221.393 -1098.342 Coefficients: Value Std.Error t-value p-value Asym.(Intercept) 44.98645 0.1958122 229.74287 0 Asym.QB300 17.24009 0.3584137 48.10108 0 lrc.(Intercept) 0.76558 0.0207001 36.98450 0 lrc.QB300 -0.51368 0.0237563 -21.62292 0 c0 0.22449 0.0028027 80.09673 0 Correlation: As.(I) A.QB30 lr.(I) l.QB30 Asym.QB300 -0.510 lrc.(Intercept) -0.649 0.256 lrc.QB300 0.513 -0.701 -0.726 c0 -0.166 -0.125 0.457 -0.081 Standardized residuals: Min Q1 Med Q3 Max -12.0842486 -1.8502594 0.3528818 2.2843122 10.2979090 Residual standard error: 1 Degrees of freedom: 140 total; 135 residual > anova (t1.fix.gnls) Denom. DF: 135 numDF F-value p-value Asym.(Intercept) 1 259134.65 <.0001 Asym.QB 1 4801.22 <.0001 lrc.(Intercept) 1 339.58 <.0001 lrc.QB 1 239.58 <.0001 c0 1 6653.10 <.0001 > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n"); .pt <- proc.time() Time elapsed: 0.02 0 0.01 > t1.fix.w <- update(t1.fix.gnls, weights = varPower()) > summary (t1.fix.w) Generalized nonlinear least squares fit Model: rate ~ SSasympOff(pressure, Asym, lrc, c0) Data: Dialyzer AIC BIC logLik 748.9766 766.6264 -368.4883 Variance function: Structure: Power of variance covariate Formula: ~fitted(.) Parameter estimates: power 0.3647475 Coefficients: Value Std.Error t-value p-value Asym.(Intercept) 45.16817 0.7723710 58.47989 0 Asym.QB300 17.52366 1.4831272 11.81535 0 lrc.(Intercept) 0.72193 0.0656526 10.99628 0 lrc.QB300 -0.50648 0.0804659 -6.29430 0 c0 0.21718 0.0047904 45.33660 0 Correlation: As.(I) A.QB30 lr.(I) l.QB30 Asym.QB300 -0.490 lrc.(Intercept) -0.668 0.271 lrc.QB300 0.501 -0.723 -0.706 c0 -0.169 -0.094 0.421 -0.083 Standardized residuals: Min Q1 Med Q3 Max -3.0472294 -0.7137222 0.1023747 0.6277747 2.6772995 Residual standard error: 1 Degrees of freedom: 140 total; 135 residual > anova (t1.fix.w) Denom. DF: 135 numDF F-value p-value Asym.(Intercept) 1 14687.148 <.0001 Asym.QB 1 249.977 <.0001 lrc.(Intercept) 1 271.105 <.0001 lrc.QB 1 6.704 0.0107 c0 1 2131.533 <.0001 > (it1fw <- intervals(t1.fix.w)) Approximate 95% confidence intervals Coefficients: lower est. upper Asym.(Intercept) 43.6406575 45.1681696 46.6956818 Asym.QB300 14.5904914 17.5236607 20.4568299 lrc.(Intercept) 0.5920936 0.7219343 0.8517749 lrc.QB300 -0.6656138 -0.5064769 -0.3473401 c0 0.2077074 0.2171813 0.2266553 Variance function: lower est. upper power 0.331479 0.3647475 0.3980161 Residual standard error: [1] 1 1 1 > stopifnot(all.equal(it1fw$varStruct["power",], + c(lower = 0.33147895, + est. = 0.36474755, + upper = 0.39801614), tol = 1e-6)) > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.05 0 0.05 > > ##=== example 4 mixed non-linear model page 363 nlme ======================= > ex <- "ex4_nlme_page363"; .pt <- proc.time() > ## > cat("\n example ", ex,"\n") example ex4_nlme_page363 > method <- "ML" > cat("\nVariable sigma; estimation method ", method,"\n") Variable sigma; estimation method ML > t1.ML.nlme <- nlme(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph, + fixed = lKe + lKa + lCl ~ 1, + method = method, start = c(-2.4,0.45,-3.2)) Warning message: In nlme.formula(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph, : Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'! > ## default control, no fixed sigma > t1.ML.nlme$numIter # 23 (32-bit) [1] 22 > stopifnot(all.equal(as.vector(t1.ML.nlme$sigma), 0.6818253, tol = 1e-5)) > > sigma <- 0.7 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 0.7 estimation method ML > set.seed(44) > system.time(# *not* fast : + t1.fix.ML.nlme <- nlme(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph, + fixed = lKe + lKa + lCl ~ 1, + method = method, start = c(-2.4,0.45,-3.2), + control = nlmeControl(sigma=sigma, maxIter = 200), + verbose = interactive()) + ) user system elapsed 3.95 0.05 4.02 Warning message: In nlme.formula(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph, : Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'! > t1.fix.ML.nlme$numIter # 58 or 61 (and now 22).. [1] 53 > (sM4 <- summary(t1.fix.ML.nlme)) Nonlinear mixed-effects model fit by maximum likelihood Model: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) Data: Theoph AIC BIC logLik 272.7945 298.7398 -127.3973 Random effects: Formula: list(lKe ~ 1, lKa ~ 1, lCl ~ 1) Level: Subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr lKe 0.1306241 lKe lKa lKa 0.6354612 0.011 lCl 0.2506970 0.995 -0.089 Residual 0.7000000 Fixed effects: lKe + lKa + lCl ~ 1 Value Std.Error DF t-value p-value lKe -2.432514 0.06401571 118 -37.99870 0.0000 lKa 0.450172 0.19605483 118 2.29615 0.0234 lCl -3.214472 0.08083845 118 -39.76415 0.0000 Correlation: lKe lKa lKa -0.149 lCl 0.852 -0.133 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.3180202 -0.4264548 0.0000000 0.2856544 3.8294100 Number of Observations: 132 Number of Groups: 12 > (aM4 <- anova (t1.fix.ML.nlme)) numDF denDF F-value p-value lKe 1 118 65.4402 <.0001 lKa 1 118 9.0874 0.0032 lCl 1 118 1581.1873 <.0001 > t1.fix.ML.nlme$apVar ## "Non-positive definite approximate variance-covariance" [1] "Non-positive definite approximate variance-covariance" > ##(iM4 <- intervals(t1.fix.ML.nlme)) > stopifnot( + all.equal(fixef(t1.fix.ML.nlme), + c(lKe = -2.432512, lKa = 0.450163, lCl = -3.2144713), tol= 8e-6) + , + all.equal(sM4$tTable[,"Std.Error"], # aarch64/linux gave 5.915e-5 + c(lKe = 0.0640155, lKa = 0.196058, lCl = 0.0808379), tol = 1e-4) + , + all.equal(aM4[,"F-value"], + c(65.439, 9.09557, 1581.21), tol = 1e-4) # ATLAS had 7.86e-05 + ) > > ## > ## REML method > if(nlme:::doExtras()) { ## -- takes 2--3 minutes + method <- "REML" + sigma <- 0.7 + cat("\nFixed sigma= ", sigma," estimation method ", method,"\n") + ## + ## only converges when tolerance is not small (and still takes long!) : + t1.fix.REML.nlme <- update(t1.fix.ML.nlme, method = method, + control = nlmeControl(tolerance = 0.0005, + sigma=sigma, + pnlsMaxIter = 20, # not just 7 + maxIter = 1000), + verbose = interactive()) + cat(" -> numIter: ", t1.fix.REML.nlme$numIter, "\n") # 380 or so + print(summary(t1.fix.REML.nlme)) + print( anova(t1.fix.REML.nlme)) + it1.fRn <- try( intervals(t1.fix.REML.nlme) ) ## cannot get .. Non-positive ... + }# only if(nlme:::doExtras()) > > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 6.45 0.08 6.55 > > ##=== example 5 mixed non-linear model page 358 nlme ======================= > ## > ## > ex <- "ex5_nlme_page365"; .pt <- proc.time() > cat("\n example ", ex,"\n") example ex5_nlme_page365 > ## > method <- "ML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method ML > t5.fix.ML.nlme <- nlme(circumference ~ SSlogis(age,Asym,xmid,scal), data = Orange, + fixed = Asym + xmid + scal ~ 1, + method = method, start = c(192,727,356), + control = nlmeControl(sigma = sigma)) > (sM5 <- summary(t5.fix.ML.nlme)) Nonlinear mixed-effects model fit by maximum likelihood Model: circumference ~ SSlogis(age, Asym, xmid, scal) Data: Orange AIC BIC logLik 1436.053 1450.051 -709.0264 Random effects: Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Level: Tree Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr Asym 30.22904 Asym xmid xmid 74.73421 0.365 scal 34.03912 -0.757 0.300 Residual 1.00000 Fixed effects: Asym + xmid + scal ~ 1 Value Std.Error DF t-value p-value Asym 192.7902 14.16878 28 13.60669 0 xmid 726.3635 35.34249 28 20.55213 0 scal 355.6294 16.36370 28 21.73283 0 Correlation: Asym xmid xmid 0.369 scal -0.722 0.315 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -11.817489 -4.184585 1.313943 4.214399 9.706200 Number of Observations: 35 Number of Groups: 5 > (aM5 <- anova (t5.fix.ML.nlme)) numDF denDF F-value p-value Asym 1 28 4879.496 <.0001 xmid 1 28 208.534 <.0001 scal 1 28 472.316 <.0001 > (t5.fix.ML.nlme$apVar) ## Non-positive definite [FIXME?] [1] "Non-positive definite approximate variance-covariance" > stopifnot( + all.equal(fixef(t5.fix.ML.nlme), + c(Asym= 192.79023, xmid= 726.36351, scal= 355.62941), tol= 1e-7) + , + all.equal(sM5$tTable[,"Std.Error"], + c(Asym= 14.1688, xmid= 35.3425, scal=16.3637), tol = 5e-5) + , + all.equal(aM5[,"F-value"], + c(4879.5, 208.534, 472.316), tol = 5e-5) + ) > > ## REML method > method <-"REML" > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method REML > t5.fix.REML.nlme <- update(t5.fix.ML.nlme, method = method, + control = nlmeControl(sigma=sigma), + verbose = interactive()) Warning message: In nlme.formula(model = circumference ~ SSlogis(age, Asym, xmid, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8) > ## converges very quickly (when started from ML!) > (sR5 <- summary(t5.fix.REML.nlme)) Nonlinear mixed-effects model fit by REML Model: circumference ~ SSlogis(age, Asym, xmid, scal) Data: Orange AIC BIC logLik 1427.774 1440.966 -704.887 Random effects: Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Level: Tree Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr Asym 30.22909 Asym xmid xmid 74.73453 0.365 scal 34.03916 -0.757 0.300 Residual 1.00000 Fixed effects: Asym + xmid + scal ~ 1 Value Std.Error DF t-value p-value Asym 192.7902 13.54797 28 14.23020 0 xmid 726.3635 33.79403 28 21.49384 0 scal 355.6294 15.64671 28 22.72870 0 Correlation: Asym xmid xmid 0.369 scal -0.722 0.315 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -11.817490 -4.184585 1.313942 4.214399 9.706200 Number of Observations: 35 Number of Groups: 5 > (aR5 <- anova (t5.fix.REML.nlme)) numDF denDF F-value p-value Asym 1 28 5336.293 <.0001 xmid 1 28 228.076 <.0001 scal 1 28 516.594 <.0001 > ( t5.fix.REML.nlme$apVar) ## Non-positive definite [FIXME?] [1] "Non-positive definite approximate variance-covariance" > stopifnot( + ## ML and REML : fixed effects are very close + all.equal(fixef(t5.fix.REML.nlme), fixef(t5.fix.ML.nlme), tol = 1e-6) + , + all.equal(sR5$tTable[,"Std.Error"], + c(Asym= 13.548, xmid= 33.794, scal=15.6467), tol = 5e-5) + , + all.equal(aR5[,"F-value"], + c(5336.29, 228.076, 516.594), tol = 5e-5) + ) > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.28 0 0.28 > > ##=== example 6 linear mixed model page 177 lme ML and REML ================ > ## > ex <- "ex6_lme_page177"; .pt <- proc.time() > cat("\n example ", ex,"\n") example ex6_lme_page177 > ## > method <- "ML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method ML > t6.fix.ML.lme <- lme(distance ~ I(age-11), data = Orthodont, + weights = varIdent(form = ~1 | Sex), + method = method, + control = lmeControl(sigma = sigma)) > (sM6 <- summary (t6.fix.ML.lme)) Linear mixed-effects model fit by maximum likelihood Data: Orthodont AIC BIC logLik 453.6113 469.7041 -220.8057 Random effects: Formula: ~I(age - 11) | Subject Structure: General positive-definite StdDev Corr (Intercept) 2.1521357 (Intr) I(age - 11) 0.2864865 0.403 Residual 1.0000000 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | Sex Parameter estimates: Male Female 1.0000000 0.6506593 Fixed effects: distance ~ I(age - 11) Value Std.Error DF t-value p-value (Intercept) 24.009565 0.4265605 80 56.28642 0 I(age - 11) 0.647604 0.0668320 80 9.69004 0 Correlation: (Intr) I(age - 11) 0.33 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.73501986 -0.57396868 0.06846467 0.53163367 5.04157906 Number of Observations: 108 Number of Groups: 27 > (aM6 <- anova (t6.fix.ML.lme)) numDF denDF F-value p-value (Intercept) 1 80 3162.474 <.0001 I(age - 11) 1 80 93.897 <.0001 > (iM6 <- intervals(t6.fix.ML.lme)) Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 23.1685797 24.0095655 24.8505512 I(age - 11) 0.5158417 0.6476043 0.7793669 Random Effects: Level: Subject lower est. upper sd((Intercept)) 1.70666362 2.1521357 2.7138847 sd(I(age - 11)) 0.20110174 0.2864865 0.4081242 cor((Intercept),I(age - 11)) 0.04335791 0.4031354 0.6703559 Variance function: lower est. upper Female 0.512246 0.6506593 0.8264729 > stopifnot( + all.equal(fixef(t6.fix.ML.lme), + c("(Intercept)"= 24.009565, "I(age - 11)"= 0.64760432), tol= 1e-7) + , + all.equal(sM6$tTable[,"Std.Error"], + c("(Intercept)"= 0.426561, "I(age - 11)"= 0.066832), tol = 5e-5) + , + all.equal(aM6[,"F-value"], c(3162.47, 93.8969), tol = 5e-5) + , + all.equal(iM6$varStruct["Female",], ## Win 32 + c(lower = 0.51230063, ## 0.51226722 + est. = 0.65065925, ## 0.65065925 + upper = 0.82638482), ## 0.82643872 + ## was tol = if(isSun) 4e-4 else 6e-5)#= 4.39e-5 + tol = 4e-4) + ## seen 5.35e-5 (Sparc Sol., no long double); later, 6e-5 was not ok + ## Windows 64bit w/ OpenBLAS 0.2.18 gave 5.721e-05 (Avi A) + ## Win-builder (i386) gave 9.584157e-05 + ) > > ##------------- > method <- "REML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method REML > > t6.fix.REML.lme <- lme(distance ~I(age-11), data = Orthodont, + weights = varIdent(form = ~1 | Sex), + method = method, + control = lmeControl(sigma = sigma)) > (sR6 <- summary (t6.fix.REML.lme)) Linear mixed-effects model fit by REML Data: Orthodont AIC BIC logLik 457.6177 473.5983 -222.8089 Random effects: Formula: ~I(age - 11) | Subject Structure: General positive-definite StdDev Corr (Intercept) 2.2247716 (Intr) I(age - 11) 0.3013982 0.394 Residual 1.0000000 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | Sex Parameter estimates: Male Female 1.000000 0.660878 Fixed effects: distance ~ I(age - 11) Value Std.Error DF t-value p-value (Intercept) 24.01066 0.4363651 80 55.02425 0 I(age - 11) 0.64880 0.0687549 80 9.43642 0 Correlation: (Intr) I(age - 11) 0.326 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.7290540 -0.5654783 0.0772640 0.5224078 5.0325191 Number of Observations: 108 Number of Groups: 27 > (aR6 <- anova (t6.fix.REML.lme)) numDF denDF F-value p-value (Intercept) 1 80 3019.859 <.0001 I(age - 11) 1 80 89.046 <.0001 > (iR6 <- intervals(t6.fix.REML.lme)) Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 23.1422680 24.0106621 24.8790563 I(age - 11) 0.5119731 0.6487997 0.7856262 Random Effects: Level: Subject lower est. upper sd((Intercept)) 1.74897367 2.2247716 2.8300075 sd(I(age - 11)) 0.21182414 0.3013982 0.4288505 cor((Intercept),I(age - 11)) 0.02154616 0.3935633 0.6698491 Variance function: lower est. upper Female 0.5177462 0.660878 0.8435788 > stopifnot( + all.equal(fixef(t6.fix.REML.lme), + c("(Intercept)"= 24.010662, "I(age - 11)"= 0.64879966), tol= 1e-7) + , + all.equal(sR6$tTable[,"Std.Error"], + c("(Intercept)"= 0.436365, "I(age - 11)"= 0.0687549), tol = 5e-5) + , + all.equal(aR6[,"F-value"], c(3019.86, 89.046), tol = 5e-5) + , + all.equal(iR6$varStruct["Female",], ## Win 32 + c(lower = 0.51774671, ## 0.51778038 + est. = 0.66087796, ## 0.66087807 + upper = 0.8435779), ## 0.84352331 + ## was tol = if(isSun) 4e-4 else 5e-5)# 4.37e-5 + tol = 4e-4) # was 5.63e-5 without long doubles. + ) > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.13 0.02 0.14 > > ##=== example 7 linear mixed model page 172 lme ML and REML ================ > ex <- "ex7_lme_page172"; .pt <- proc.time() > > cat("\n example ", ex,"\n") example ex7_lme_page172 > method <- "ML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method ML > set.seed(107) > t7.fix.ML.lme <- lme( current ~ voltage + I(voltage^2), data = Wafer, + random = list(Wafer = pdDiag(~voltage + I(voltage^2)), + Site = pdDiag(~voltage + I(voltage^2)) ), + method = method, + control = lmeControl(sigma = 1, + ## nlminb: false convergence on 32-bit + msVerbose = TRUE, opt = "optim")) initial value 30.453061 iter 10 value 23.916897 final value 23.913475 converged > (ss7 <- summary(t7.fix.ML.lme)) Linear mixed-effects model fit by maximum likelihood Data: Wafer AIC BIC logLik 800.9778 836.901 -391.4889 Random effects: Formula: ~voltage + I(voltage^2) | Wafer Structure: Diagonal (Intercept) voltage I(voltage^2) StdDev: 0.001097033 0.1620467 0.0008636442 Formula: ~voltage + I(voltage^2) | Site %in% Wafer Structure: Diagonal (Intercept) voltage I(voltage^2) Residual StdDev: 0.002149923 0.001139061 0.0005960515 1 Fixed effects: current ~ voltage + I(voltage^2) Value Std.Error DF t-value p-value (Intercept) -4.461166 0.4460859 318 -10.000688 0 voltage 5.903371 0.6085713 318 9.700377 0 I(voltage^2) 1.170403 0.1874590 318 6.243514 0 Correlation: (Intr) voltag voltage -0.974 I(voltage^2) 0.941 -0.986 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -0.62903588 -0.18062130 -0.01551271 0.13500732 1.05129039 Number of Observations: 400 Number of Groups: Wafer Site %in% Wafer 10 80 > (aa7 <- anova(t7.fix.ML.lme)) numDF denDF F-value p-value (Intercept) 1 318 2634.214 <.0001 voltage 1 318 8851.051 <.0001 I(voltage^2) 1 318 38.981 <.0001 > stopifnot( + all.equal(fixef(t7.fix.ML.lme), + c("(Intercept)" = -4.4611657, "voltage" = 5.9033709, + "I(voltage^2)" = 1.1704027), tol = 1e-7) + , + all.equal(ss7$tTable[,"Std.Error"], + c("(Intercept)" = 0.446086, "voltage" = 0.608571, + "I(voltage^2)"= 0.187459), tol = 5e-5) + , + all.equal(aa7[,"F-value"], + c(2634.2137, 8851.0513, 38.981463), tol = 5e-5) + ) > ##------------------------------------------------------ REML --- > method <- "REML" > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method REML > ## This had 'false convergence' on Solaris using nlminb, and > ## optim found a portable answer > t7.fix.REML.lme <- lme( current ~ voltage + I(voltage^2), data = Wafer, + random = list(Wafer = pdDiag(~voltage + I(voltage^2)), + Site = pdDiag(~voltage + I(voltage^2)) ), + control = lmeControl(sigma = 1, opt = "optim"), + method = method) > (sR7 <- summary(t7.fix.REML.lme)) Linear mixed-effects model fit by REML Data: Wafer AIC BIC logLik 801.5136 837.3691 -391.7568 Random effects: Formula: ~voltage + I(voltage^2) | Wafer Structure: Diagonal (Intercept) voltage I(voltage^2) StdDev: 0.001444446 0.1632295 0.0009738082 Formula: ~voltage + I(voltage^2) | Site %in% Wafer Structure: Diagonal (Intercept) voltage I(voltage^2) Residual StdDev: 0.00235956 0.001229668 0.0006564859 1 Fixed effects: current ~ voltage + I(voltage^2) Value Std.Error DF t-value p-value (Intercept) -4.461166 0.4444100 318 -10.038400 0 voltage 5.903371 0.6063166 318 9.736450 0 I(voltage^2) 1.170403 0.1867547 318 6.267057 0 Correlation: (Intr) voltag voltage -0.974 I(voltage^2) 0.941 -0.986 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -0.62763794 -0.18117128 -0.01568026 0.13501724 1.04998361 Number of Observations: 400 Number of Groups: Wafer Site %in% Wafer 10 80 > (aR7 <- anova(t7.fix.REML.lme)) numDF denDF F-value p-value (Intercept) 1 318 2584.952 <.0001 voltage 1 318 8885.000 <.0001 I(voltage^2) 1 318 39.276 <.0001 > stopifnot( + all.equal(fixef(t7.fix.REML.lme), fixef(t7.fix.ML.lme), + ## should not change much from ML to REML ! + tol = 1e-6) + , + all.equal(sR7$tTable[,"Std.Error"], + c("(Intercept)" = 0.44441, "voltage" = 0.606321, + "I(voltage^2)"= 0.186754), tol = 5e-5) + , + all.equal(aR7[,"F-value"], + c(2584.9515, 8885.0, 39.2760), tol = 1e-6) + ) > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 1.62 0 1.62 > > ##=== example 8 mixed non-linear model page 364 nlme ======================= > ## > ex <- "ex8_nlme_page364"; .pt <- proc.time() > ## > cat("\n example ", ex,"\n") example ex8_nlme_page364 > method <- "ML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method ML > set.seed(8^2) > t8.fix.ML.nlme <- nlme(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph, + fixed = lKe + lKa + lCl ~ 1, + random = pdDiag(lKe + lKa + lCl ~ 1), + method = method, start = c(-2.4,0.5,-3.3), + control = nlmeControl(sigma = 1)) > (sM8 <- summary(t8.fix.ML.nlme)) Nonlinear mixed-effects model fit by maximum likelihood Model: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) Data: Theoph AIC BIC logLik 386.8467 404.1435 -187.4234 Random effects: Formula: list(lKe ~ 1, lKa ~ 1, lCl ~ 1) Level: Subject Structure: Diagonal lKe lKa lCl Residual StdDev: 6.424297e-06 0.6072879 0.1621444 1 Fixed effects: lKe + lKa + lCl ~ 1 Value Std.Error DF t-value p-value lKe -2.455500 0.07392689 118 -33.21525 0.0000 lKa 0.448703 0.19752391 118 2.27164 0.0249 lCl -3.229696 0.06830493 118 -47.28349 0.0000 Correlation: lKe lKa lKa -0.270 lCl 0.668 -0.135 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.2873732 -0.3872350 0.0000000 0.1946299 2.7147893 Number of Observations: 132 Number of Groups: 12 > (aM8 <- anova (t8.fix.ML.nlme)) numDF denDF F-value p-value lKe 1 118 10.9426 0.0012 lKa 1 118 17.4101 0.0001 lCl 1 118 2235.7286 <.0001 > stopifnot( + all.equal(fixef(t8.fix.ML.nlme), + c(lKe = -2.4554999, lKa = 0.44870292, lCl = -3.2296957), + tol = 1e-7) + , + all.equal(sM8$tTable[,"Std.Error"], + c(lKe = 0.0739269, lKa = 0.197524, lCl = 0.0683049), + tol = 5e-5) + , + all.equal(aM8[,"F-value"], + c(10.9426, 17.4101, 2235.73), tol = 5e-5) + ) > ## REML method > method <- "REML" > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method REML > t8.fix.REML.nlme <- update(t8.fix.ML.nlme, method = method) > (sR8 <- summary(t8.fix.REML.nlme)) Nonlinear mixed-effects model fit by REML Model: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) Data: Theoph AIC BIC logLik 389.6652 406.8241 -188.8326 Random effects: Formula: list(lKe ~ 1, lKa ~ 1, lCl ~ 1) Level: Subject Structure: Diagonal lKe lKa lCl Residual StdDev: 8.696678e-06 0.6072879 0.1621444 1 Fixed effects: lKe + lKa + lCl ~ 1 Value Std.Error DF t-value p-value lKe -2.455500 0.07308198 118 -33.59925 0.0000 lKa 0.448703 0.19526642 118 2.29790 0.0233 lCl -3.229696 0.06752428 118 -47.83014 0.0000 Correlation: lKe lKa lKa -0.270 lCl 0.668 -0.135 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.2873732 -0.3872350 0.0000000 0.1946299 2.7147893 Number of Observations: 132 Number of Groups: 12 > (aR8 <- anova (t8.fix.REML.nlme)) numDF denDF F-value p-value lKe 1 118 11.1971 0.0011 lKa 1 118 17.8150 <.0001 lCl 1 118 2287.7223 <.0001 > stopifnot( + all.equal(fixef(t8.fix.REML.nlme), fixef(t8.fix.ML.nlme), + ## should not change much from ML to REML ! + tol = 1e-6) + , + all.equal(sR8$tTable[,"Std.Error"], + c(lKe = 0.073082, lKa = 0.195266, lCl = 0.0675243), + tol = 5e-5) + , + all.equal(aR8[,"F-value"], + c(11.1971, 17.815, 2287.72), tol = 5e-5) + ) > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.24 0.01 0.25 > > ##=== example 9 mixed non-linear model page 365 nlme ======================= > ## > ex <- "ex9_nlme_page365"; .pt <- proc.time() > ## > cat("\n example ", ex,"\n") example ex9_nlme_page365 > method <- "ML" > sigma <- 1 > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method ML > set.seed(909) > t9.fix.ML.nlme <- nlme(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph, + fixed = lKe + lKa + lCl ~ 1, + random = pdDiag( lKa + lCl ~ 1), + method = method, start = c(-2.4,0.5,-3.3), + control = nlmeControl(sigma = 1)) > (sM9 <- summary(t9.fix.ML.nlme)) Nonlinear mixed-effects model fit by maximum likelihood Model: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) Data: Theoph AIC BIC logLik 384.8436 399.2576 -187.4218 Random effects: Formula: list(lKa ~ 1, lCl ~ 1) Level: Subject Structure: Diagonal lKa lCl Residual StdDev: 0.6070128 0.1621536 1 Fixed effects: lKe + lKa + lCl ~ 1 Value Std.Error DF t-value p-value lKe -2.455574 0.07392660 118 -33.21639 0.0000 lKa 0.448941 0.19745878 118 2.27359 0.0248 lCl -3.229727 0.06830818 118 -47.28171 0.0000 Correlation: lKe lKa lKa -0.270 lCl 0.668 -0.136 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.2877079 -0.3872728 0.0000000 0.1947776 2.7147395 Number of Observations: 132 Number of Groups: 12 > (aM9 <- anova (t9.fix.ML.nlme)) numDF denDF F-value p-value lKe 1 118 10.9669 0.0012 lKa 1 118 17.4108 0.0001 lCl 1 118 2235.5599 <.0001 > stopifnot( + all.equal(fixef(t9.fix.ML.nlme), + c(lKe = -2.4555745, lKa = 0.44894103, lCl = -3.2297273), + tol = 1e-7) + , + all.equal(sM9$tTable[,"Std.Error"], + c(lKe = 0.0739266, lKa = 0.197459, lCl = 0.0683082), + tol = 5e-5) + , + all.equal(aM9[,"F-value"], + c(10.9669, 17.4108, 2235.56), tol = 5e-5) + ) > ## REML method > method <- "REML" > cat("\nFixed sigma= ",sigma," estimation method ", method,"\n") Fixed sigma= 1 estimation method REML > t9.fix.REML.nlme <- update(t9.fix.ML.nlme, method = method) > (sR9 <- summary(t9.fix.REML.nlme)) Nonlinear mixed-effects model fit by REML Model: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) Data: Theoph AIC BIC logLik 387.6629 401.962 -188.8315 Random effects: Formula: list(lKa ~ 1, lCl ~ 1) Level: Subject Structure: Diagonal lKa lCl Residual StdDev: 0.6070128 0.1621536 1 Fixed effects: lKe + lKa + lCl ~ 1 Value Std.Error DF t-value p-value lKe -2.455574 0.07308170 118 -33.60040 0.0000 lKa 0.448941 0.19520204 118 2.29988 0.0232 lCl -3.229727 0.06752749 118 -47.82834 0.0000 Correlation: lKe lKa lKa -0.270 lCl 0.668 -0.136 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.2877079 -0.3872728 0.0000000 0.1947776 2.7147395 Number of Observations: 132 Number of Groups: 12 > (aR9 <- anova (t9.fix.REML.nlme)) numDF denDF F-value p-value lKe 1 118 11.2219 0.0011 lKa 1 118 17.8157 <.0001 lCl 1 118 2287.5497 <.0001 > stopifnot( + all.equal(fixef(t9.fix.REML.nlme), fixef(t9.fix.ML.nlme), + ## should not change much from ML to REML ! + tol = 1e-6) + , + all.equal(sR9$tTable[,"Std.Error"], + c(lKe = 0.0730817, lKa = 0.195202, lCl = 0.0675275), + tol = 5e-5) + , + all.equal(aR9[,"F-value"], + c(11.2219, 17.8157, 2287.55), tol = 5e-5) + ) > cat("Time elapsed: ", (proc.time() - .pt)[1:3], "\n") Time elapsed: 0.14 0 0.14 > > proc.time() user system elapsed 9.53 0.18 9.71