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("segmented-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('segmented')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "down"
> 
> ### * down
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: down
> ### Title: Down syndrome in babies
> ### Aliases: down
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(down)
> fit.glm<-glm(cases/births~age, weight=births, family=binomial, data=down)
> fit.seg<-segmented(fit.glm, Z=age, psi=25)
> fit.seg
Call: segmented.glm(obj = fit.glm, Z = age, psi = 25)

Meaningful coefficients of the linear terms:
(Intercept)         age       U.age 
-6.78243778 -0.01341037  0.27422124 

Estimated Break-Point(s) for the variable(s) age : 31.08 

Degrees of Freedom: 29 Total (i.e. Null);  26 Residual
Null Deviance:     625.2 
Residual Deviance: 43.94      AIC: 190.8 
> 
> 
> 
> cleanEx(); ..nameEx <- "segmented"
> 
> ### * segmented
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: segmented
> ### Title: Segmented relationships in regression models
> ### Aliases: segmented segmented.default segmented.lm segmented.glm
> ###   print.segmented summary.segmented print.summary.segmented
> ### Keywords: regression
> 
> ### ** Examples
> 
> ## A linear model with a segmented relationship in each level of 
> ## a categorical variable
> ##
> ## Simulate data
> x<-1:100
> g<-rep(0:1,c(100,100))
> set.seed(15)
> y1<- 2+1.5*pmax(x-60,0)+rnorm(100,0,5)
> y2<- 2+0.6*pmax(x-60,0)+rnorm(100,0,5)
> dati<-data.frame(xx=rep(x,2),yy=c(y1,y2),g=factor(g))
> rm(x,g,y1,y2)
> ## Have a look at the plot
> ## Not run: plot(dati$xx,dati$yy) 
> ## Not run: points(dati$xx[dati$g==0],dati$yy[dati$g==0],col=2)
> ## Fit the linear model
> obj<-lm(yy~0+g+xx:g,data=dati) 
> ## Fit segmented models
> ## Model I: ignore the stratification, assume an equal difference-in-slopes 
> ## parameter between groups
> ogg0<-segmented(obj,Z=xx,psi=60)
> ## Model II: now stratificate by g. Here note that psi[i] refers to the 
> ## segmented relationship in the ith level of W
> ogg1<-segmented(obj,Z=xx,W=g,psi=c(50,70),visual=TRUE)
0   5993.535  (No breakpoint(s)) 
1   5993.535 
2   5223.672 
3   5146.034 
4   5146.034 
> ## Results....
> ogg0
Call: segmented.lm(obj = obj, Z = xx, psi = 60)

Meaningful coefficients of the linear terms:
           g0            g1          U.xx        psi.xx         g0:xx 
-4.974130e-01  5.893234e+00  1.119724e+00  9.916040e-16  1.285178e-01 

Estimated Break-Point(s) for the variable(s) xx : 60.81 
> summary(ogg1)

	***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = obj, Z = xx, psi = c(50, 70), W = g, visual = TRUE)

Estimated Break-Point(s):
        Est. St.Err
g0:xx 63.24  1.362
g1:xx 57.44  2.976

t value for the gap-variable(s) V:  0.633222 -1.32283 

Meaningful coefficients of the linear terms:
                  Estimate Std. Error       t value
g0            2.689815e+00 1.32019009  2.037445e+00
g1            2.623828e+00 1.38969088  1.888066e+00
`U1.g0:xx`    1.581413e+00 0.08741195  1.809149e+01
`U2.g1:xx`    7.032423e-01 0.07605725  9.246223e+00
`psi1.g0:xx`  1.241564e-15 1.36227841  9.113881e-16
`psi2.g1:xx` -6.521647e-15 2.97589713 -2.191489e-15


Residual standard error: 5.177 on 192 degrees of freedom
Multiple R-Squared: 0.9272,  Adjusted R-squared: 0.9242 

Convergence attained in 4 iterations with relative change 1.76737e-16 
> ## Have a look at the fitted models
> ## model I
> ## Not run: points(dati$xx[dati$g==0],ogg0$fitted[dati$g==0],pch=3,col=2)
> ## Not run: points(dati$xx[dati$g==1],ogg0$fitted[dati$g==1],pch=3,col=1)
> ## model II
> ## Not run: points(dati$xx[dati$g==0],ogg1$fitted[dati$g==0],pch=20,col=2)
> ## Not run: points(dati$xx[dati$g==1],ogg1$fitted[dati$g==1],pch=20,col=1)
> ## Not run: legend(5,60,legend=c("model I","model II"),pch=c(3,20))
> 
> 
> 
> cleanEx(); ..nameEx <- "slope.segmented"
> 
> ### * slope.segmented
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: slope.segmented
> ### Title: Summary for slopes of segmented relationships
> ### Aliases: slope.segmented
> ### Keywords: regression
> 
> ### ** Examples
> 
> ## A segmented relationship with three breakpoints
> ## 
> x<-1:100
> y<-2+1.5*pmax(x-35,0)-1.5*pmax(x-70,0)+rnorm(100,0,5)
> dati<-data.frame(x=x,y=y)
> out.lm<-lm(y~x,data=dati)
> out.seg<-segmented(out.lm,Z=x,psi=c(20,80))
Warning: max number of iterations attained
> ## the slopes of the three segments....
> slope.segmented(out.seg)
$x
              Est.    St.Err.    t value  CI(95%).l CI(95%).u
slope1 -0.04845980 0.07626608 -0.6354043 -0.1979386 0.1010190
slope2  1.47854910 0.08724557 16.9469824  1.3075509 1.6495473
slope3 -0.03122188 0.08330763 -0.3747782 -0.1945018 0.1320581

> rm(x,y,out.lm,out.seg)
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  0.61 0.04 0.65 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')