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("perturb-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('perturb')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "colldiag"
> 
> ### * colldiag
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: colldiag
> ### Title: Condition indexes and variance decomposition proportions
> ### Aliases: colldiag print.colldiag
> ### Keywords: regression
> 
> ### ** Examples
> 
> # Belsley (1991). "Conditioning Diagnostics"
> # The Consumption Function (pp. 149-154)
> data(consumption)
> attach(consumption)
> 
> ct1<-c(NA,c[-length(c)])
> 
> # compare (5.3)
> m1<-lm(c~ct1+dpi+r+d_dpi)
> anova(m1)
Analysis of Variance Table

Response: c
          Df Sum Sq Mean Sq    F value    Pr(>F)    
ct1        1 300165  300165 23721.4022 < 2.2e-16 ***
dpi        1   1653    1653   130.6481 1.006e-10 ***
r          1     24      24     1.9073    0.1811    
d_dpi      1     10      10     0.7685    0.3902    
Residuals 22    278      13                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary(m1)

Call:
lm(formula = c ~ ct1 + dpi + r + d_dpi)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5567 -2.5185 -0.8726  2.2804  5.8832 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   6.7242     3.8271   1.757  0.09283 . 
ct1           0.2454     0.2375   1.033  0.31271   
dpi           0.6984     0.2077   3.363  0.00281 **
r            -2.2097     1.8384  -1.202  0.24217   
d_dpi         0.1608     0.1834   0.877  0.39016   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 3.557 on 22 degrees of freedom
Multiple R-Squared: 0.9991,	Adjusted R-squared: 0.9989 
F-statistic:  5964 on 4 and 22 DF,  p-value: < 2.2e-16 

> 
> # compare exhibit 5.11
> cor(cbind(ct1,dpi,r,d_dpi),use="complete.obs")
            ct1       dpi         r     d_dpi
ct1   1.0000000 0.9973816 0.9746852 0.3138422
dpi   0.9973816 1.0000000 0.9669347 0.3765246
r     0.9746852 0.9669347 1.0000000 0.2290986
d_dpi 0.3138422 0.3765246 0.2290986 1.0000000
> 
> # compare exhibit 5.12
> cd<-colldiag(m1)
> cd
Condition
Index	Variance Decomposition Proportions
           intercept ct1   dpi   r     d_dpi
1    1.000 0.201     0.212 0.213 0.208 0.167
2    4.143 0.032     0.041 0.028 0.071 0.828
3    7.799 0.742     0.033 0.043 0.182 0.001
4   39.406 0.025     0.213 0.219 0.540 0.003
5  375.614 0.000     0.501 0.498 0.000 0.001
> print(cd,fuzz=.3)
Condition
Index	Variance Decomposition Proportions
           intercept ct1   dpi   r     d_dpi
1    1.000  .         .     .     .     .   
2    4.143  .         .     .     .    0.828
3    7.799 0.742      .     .     .     .   
4   39.406  .         .     .    0.540  .   
5  375.614  .        0.501 0.498  .     .   
> 
> 
> 
> cleanEx(); ..nameEx <- "consumption"
> 
> ### * consumption
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: consumption
> ### Title: The Consumption Function
> ### Aliases: consumption
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(consumption)
> 
> 
> 
> cleanEx(); ..nameEx <- "perturb"
> 
> ### * perturb
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: perturb
> ### Title: Perturbation analysis to evaluate collinearity
> ### Aliases: perturb summary.perturb print.summary.perturb
> ### Keywords: regression
> 
> ### ** Examples
> 
> require(car)
Loading required package: car
[1] TRUE
> data(Duncan)
> attach(Duncan)
> m1<-lm(prestige~income+education)
> summary(m1)

Call:
lm(formula = prestige ~ income + education)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.5380  -6.4174   0.6546   6.6051  34.6412 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.06466    4.27194  -1.420    0.163    
income       0.59873    0.11967   5.003 1.05e-05 ***
education    0.54583    0.09825   5.555 1.73e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 13.37 on 42 degrees of freedom
Multiple R-Squared: 0.8282,	Adjusted R-squared:  0.82 
F-statistic: 101.2 on 2 and 42 DF,  p-value: < 2.2e-16 

> anova(m1)
Analysis of Variance Table

Response: prestige
          Df  Sum Sq Mean Sq F value    Pr(>F)    
income     1 30664.8 30664.8 171.570 < 2.2e-16 ***
education  1  5516.1  5516.1  30.863 1.727e-06 ***
Residuals 42  7506.7   178.7                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> vif(m1)
   income education 
 2.104900  2.104900 
> p1<-perturb(m1,pvars=c("income","education"),prange=c(1,1))
> summary(p1)
Perturb variables:
income 		 normal(0,1) 
education 		 normal(0,1) 

Impact of perturbations on coefficients:
              mean  s.d.    min    max
(Intercept) -6.011 0.312 -6.721 -5.280
income       0.597 0.010  0.566  0.634
education    0.547 0.008  0.518  0.574
> m2<-lm(prestige~income+education+type)
> summary(m2)

Call:
lm(formula = prestige ~ income + education + type)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.890  -5.740  -1.754   5.442  28.972 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.18503    3.71377  -0.050  0.96051    
income        0.59755    0.08936   6.687 5.12e-08 ***
education     0.34532    0.11361   3.040  0.00416 ** 
typeprof     16.65751    6.99301   2.382  0.02206 *  
typewc      -14.66113    6.10877  -2.400  0.02114 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 9.744 on 40 degrees of freedom
Multiple R-Squared: 0.9131,	Adjusted R-squared: 0.9044 
F-statistic:   105 on 4 and 40 DF,  p-value: < 2.2e-16 

> anova(m2)
Analysis of Variance Table

Response: prestige
          Df  Sum Sq Mean Sq F value    Pr(>F)    
income     1 30664.8 30664.8 322.962 < 2.2e-16 ***
education  1  5516.1  5516.1  58.096 2.590e-09 ***
type       2  3708.7  1854.4  19.530 1.208e-06 ***
Residuals 40  3798.0    94.9                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> vif(m2)
              GVIF Df GVIF^(1/2Df)
income    2.209178  1     1.486330
education 5.297584  1     2.301648
type      5.098592  2     1.502666
> p2<-perturb(m2,pvars=c("income","education"),prange=c(1,1),pfac=list("type",pcnt=95))
> summary(p2)
Perturb variables:
income 		 normal(0,1) 
education 		 normal(0,1) 


Variable type to be reclassified Reclassification probabilities:
original	reclassified
        bc  prof    wc
bc   0.963 0.023 0.013
prof 0.027 0.958 0.014
wc   0.046 0.043 0.911

Impact of perturbations on coefficients:
               mean  s.d.     min    max
(Intercept)  -1.825 2.004  -7.962  2.486
income        0.602 0.042   0.423  0.721
education     0.409 0.074   0.232  0.595
typeprof     12.257 5.533  -6.036 24.326
typewc      -15.430 4.085 -23.915 -4.076
> 
> ## Not run: 
> ##D r1<-lm(ses~fegp6+educyr+eyr+exp2)
> ##D summary(r1)
> ##D q<-perturb(r1,c("eyr","exp"),c(2.5,2.5),ptrans="exp2<-exp^2")
> ##D summary(q)
> ##D 
> ##D fegp6<-as.factor(fegp6)
> ##D 
> ##D # eyr and exp also as factors
> ##D eyrc<-cut(eyr,c(min(eyr),40,50,60,70,80,max(eyr)),include.lowest=T,right=F)
> ##D table(eyrc)
> ##D expc<-cut(exp,c(0,10,20,max(exp)),include.lowest=T,right=F)
> ##D table(expc)
> ##D 
> ##D # rough initial reclassification probabilities,
> ##D # program will ensure they sum to 100 row-wise
> ##D m1<-matrix(0,nlevels(eyrc),nlevels(eyrc))
> ##D m1[row(m1)==col(m1)]<-80
> ##D m1[abs(row(m1)-col(m1))==1]<-8
> ##D m1[abs(row(m1)-col(m1))==2]<-2
> ##D m1
> ##D 
> ##D m2<-matrix(0,nlevels(expc),nlevels(expc))
> ##D m2[row(m2)==col(m2)]<-80
> ##D m2[abs(row(m2)-col(m2))==1]<-10
> ##D m2[abs(row(m2)-col(m2))==2]<-2
> ##D m2
> ##D 
> ##D r2<-lm(ses~fegp6+eyrc+expc)
> ##D summary(r2)
> ##D pfac=list(list("fegp6",pcnt=95),list("eyrc",pcnt=m1),list("expc",pcnt=m2))
> ##D q2<-perturb(r2,pfac=pfac,niter=1)
> ##D summary(q2)
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "reclassify"
> 
> ### * reclassify
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: reclassify
> ### Title: Called by perturb to calculate reclassification tables
> ### Aliases: reclassify print.reclassify
> ### Keywords: regression
> 
> ### ** Examples
> 
> library(car)
> data(Duncan)
> attach(Duncan)
> 
> reclassify(type,pcnt=95)

Variable type to be reclassified Reclassification probabilities:
original	reclassified
        bc  prof    wc
bc   0.963 0.023 0.013
prof 0.027 0.958 0.014
wc   0.046 0.043 0.911

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