R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

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 a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("lmeSplines-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('lmeSplines')
Loading required package: nlme
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "approx.Z"
> 
> ### * approx.Z
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: approx.Z
> ### Title: Interpolating in smoothing spline Z-matrix columns
> ### Aliases: approx.Z
> ### Keywords: manip
> 
> ### ** Examples
> 
> times1 <- 1:10
> Zt1 <- smspline(~ times1)
> times2 <- seq(1,10,by=0.1)
> Zt2 <- approx.Z(Zt1,oldtimes=times1,newtimes=times2)
> 
> 
> 
> cleanEx(); ..nameEx <- "smSplineEx1"
> 
> ### * smSplineEx1
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: smSplineEx1
> ### Title: Simulated data about a smooth curve
> ### Aliases: smSplineEx1
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(smSplineEx1)
> 
> 
> 
> cleanEx(); ..nameEx <- "smspline"
> 
> ### * smspline
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: smspline
> ### Title: Smoothing splines in NLME
> ### Aliases: smspline smspline.v
> ### Keywords: smooth models regression
> 
> ### ** Examples
> 
> # smoothing spline curve fit
> data(smSplineEx1)
> # variable `all' for top level grouping
> smSplineEx1$all <- rep(1,nrow(smSplineEx1))
> # setup spline Z-matrix
> smSplineEx1$Zt <- smspline(~ time, data=smSplineEx1)
> fit1s <- lme(y ~ time, data=smSplineEx1,
+     random=list(all=pdIdent(~Zt - 1)))
> summary(fit1s)
Linear mixed-effects model fit by REML
 Data: smSplineEx1 
       AIC      BIC    logLik
  281.9241 292.2640 -136.9620

Random effects:
 Formula: ~Zt - 1 | all
 Structure: Multiple of an Identity
               Zt1        Zt2        Zt3        Zt4        Zt5        Zt6
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
               Zt7        Zt8        Zt9       Zt10       Zt11       Zt12
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt13       Zt14       Zt15       Zt16       Zt17       Zt18
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt19       Zt20       Zt21       Zt22       Zt23       Zt24
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt25       Zt26       Zt27       Zt28       Zt29       Zt30
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt31       Zt32       Zt33       Zt34       Zt35       Zt36
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt37       Zt38       Zt39       Zt40       Zt41       Zt42
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt43       Zt44       Zt45       Zt46       Zt47       Zt48
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt49       Zt50       Zt51       Zt52       Zt53       Zt54
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt55       Zt56       Zt57       Zt58       Zt59       Zt60
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt61       Zt62       Zt63       Zt64       Zt65       Zt66
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt67       Zt68       Zt69       Zt70       Zt71       Zt72
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt73       Zt74       Zt75       Zt76       Zt77       Zt78
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt79       Zt80       Zt81       Zt82       Zt83       Zt84
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt85       Zt86       Zt87       Zt88       Zt89       Zt90
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt91       Zt92       Zt93       Zt94       Zt95       Zt96
StdDev: 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551 0.01627551
              Zt97       Zt98  Residual
StdDev: 0.01627551 0.01627551 0.8603034

Fixed effects: y ~ time 
               Value  Std.Error DF  t-value p-value
(Intercept) 6.498800 0.17335926 98 37.48747       0
time        0.038723 0.00298033 98 12.99299       0
 Correlation: 
     (Intr)
time -0.868

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.8893070 -0.5811366  0.1164528  0.6593223  1.7839013 

Number of Observations: 100
Number of Groups: 1 
> plot(smSplineEx1$time,smSplineEx1$y,pch="o",type="n",
+      main="Spline fits: lme(y ~ time, random=list(all=pdIdent(~Zt-1)))",
+      xlab="time",ylab="y")
> points(smSplineEx1$time,smSplineEx1$y,col=1)
> lines(smSplineEx1$time, smSplineEx1$y.true,col=1)
> lines(smSplineEx1$time, fitted(fit1s),col=2)
> 
> # fit model with cut down number of spline points
> times20 <- seq(1,100,length=20)
> Zt20 <- smspline(times20)
> smSplineEx1$Zt20 <- approx.Z(Zt20,times20,smSplineEx1$time)
> fit1s20 <- lme(y ~ time, data=smSplineEx1,
+     random=list(all=pdIdent(~Zt20 - 1)))
> # note: virtually identical df, loglik.
> anova(fit1s,fit1s20)
        Model df      AIC      BIC    logLik
fit1s       1  4 281.9241 292.2640 -136.9621
fit1s20     2  4 281.9491 292.2889 -136.9745
> summary(fit1s20)
Linear mixed-effects model fit by REML
 Data: smSplineEx1 
       AIC      BIC    logLik
  281.9491 292.2889 -136.9745

Random effects:
 Formula: ~Zt20 - 1 | all
 Structure: Multiple of an Identity
             Zt201      Zt202      Zt203      Zt204      Zt205      Zt206
StdDev: 0.01613012 0.01613012 0.01613012 0.01613012 0.01613012 0.01613012
             Zt207      Zt208      Zt209     Zt2010     Zt2011     Zt2012
StdDev: 0.01613012 0.01613012 0.01613012 0.01613012 0.01613012 0.01613012
            Zt2013     Zt2014     Zt2015     Zt2016     Zt2017     Zt2018
StdDev: 0.01613012 0.01613012 0.01613012 0.01613012 0.01613012 0.01613012
         Residual
StdDev: 0.8618987

Fixed effects: y ~ time 
               Value  Std.Error DF  t-value p-value
(Intercept) 6.392423 0.17557145 98 36.40924       0
time        0.039741 0.00302297 98 13.14633       0
 Correlation: 
     (Intr)
time -0.87 

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.8961975 -0.5735940  0.1249486  0.6647808  1.7871069 

Number of Observations: 100
Number of Groups: 1 
> 
> # model predictions on a finer grid
> times200 <- seq(1,100,by=0.5)
> pred.df <- data.frame(all=rep(1,length(times200)),time=times200)
> pred.df$Zt20 <- approx.Z(Zt20, times20,times200)
> yp20.200 <- predict(fit1s20,newdata=pred.df)
> lines(times200,yp20.200+0.02,col=4)
> 
> # mixed model spline terms at multiple levels of grouping
> data(Spruce)
> Spruce$Zday <- smspline(~ days, data=Spruce)
> Spruce$all <- rep(1,nrow(Spruce))
> # overall spline term, random plot and Tree effects
> spruce.fit1 <- lme(logSize ~ days, data=Spruce,
+                    random=list(all= pdIdent(~Zday -1),
+                      plot=~1, Tree=~1))
> # try overall spline term plus plot level linear + spline term
> spruce.fit2 <- lme(logSize ~ days, data=Spruce,
+                    random=list(all= pdIdent(~Zday - 1),
+                      plot= pdBlocked(list(~ days,pdIdent(~Zday - 1))),
+                      Tree = ~1))
Warning: Singular precision matrix in level -2, block 1
Warning: Singular precision matrix in level -2, block 1
Warning: Singular precision matrix in level -2, block 1
> anova(spruce.fit1,spruce.fit2)
            Model df       AIC       BIC    logLik   Test  L.Ratio p-value
spruce.fit1     1  6 -178.6251 -149.0304  95.31254                        
spruce.fit2     2  9 -230.3817 -185.9897 124.19085 1 vs 2 57.75661  <.0001
> summary(spruce.fit1)
Linear mixed-effects model fit by REML
 Data: Spruce 
        AIC       BIC   logLik
  -178.6251 -149.0304 95.31254

Random effects:
 Formula: ~Zday - 1 | all
 Structure: Multiple of an Identity
               Zday1        Zday2        Zday3        Zday4        Zday5
StdDev: 0.0007578229 0.0007578229 0.0007578229 0.0007578229 0.0007578229
               Zday6        Zday7        Zday8        Zday9       Zday10
StdDev: 0.0007578229 0.0007578229 0.0007578229 0.0007578229 0.0007578229
              Zday11
StdDev: 0.0007578229

 Formula: ~1 | plot %in% all
        (Intercept)
StdDev:  0.08899166

 Formula: ~1 | Tree %in% plot %in% all
        (Intercept)  Residual
StdDev:   0.6217336 0.1765053

Fixed effects: logSize ~ days 
               Value  Std.Error  DF   t-value p-value
(Intercept) 4.140108 0.08520752 947  48.58853       0
days        0.003322 0.00002940 947 112.99156       0
 Correlation: 
     (Intr)
days -0.148

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-4.433904446 -0.482001090 -0.007196469  0.535761025  6.674696431 

Number of Observations: 1027
Number of Groups: 
                    all           plot %in% all Tree %in% plot %in% all 
                      1                       4                      79 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  5.37 0.51 5.94 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')