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("psy-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('psy')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "ckappa"
> 
> ### * ckappa
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ckappa
> ### Title: Cohen's Kappa for 2 raters
> ### Aliases: ckappa
> ### Keywords: univar
> 
> ### ** Examples
> data(expsy)
> ckappa(expsy[,c(12,14)])          # Cohen's kappa for binary diagnosis
$table
  0  1
0 6  0
1 4 19

$kappa
[1] 0.6627907

> 
> library(boot)
> ckappa.boot <- function(data,x) {ckappa(data[x,])[[2]]}
> res <- boot(expsy[,c(12,14)],ckappa.boot,1000)
> quantile(res$t,c(0.025,0.975))    # two-sided bootstrapped confidence interval of kappa
     2.5%     97.5% 
0.3547560 0.9204890 
> boot.ci(res,type="bca")           # adjusted bootstrap percentile (BCa) confidence interval (better)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = "bca")

Intervals : 
Level       BCa          
95%   ( 0.3100,  0.9112 )  
Calculations and Intervals on Original Scale
> ckappa(expsy[,c(11,13)])          # Cohen's kappa for non binary diagnosis
$table
  1 2 3 4
1 1 0 0 0
2 1 4 0 0
3 0 4 9 2
4 0 0 2 6

$kappa
[1] 0.5421053

> 
> 
> 
> cleanEx(); ..nameEx <- "cronbach"
> 
> ### * cronbach
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cronbach
> ### Title: Cronbach's coefficient alpha
> ### Aliases: cronbach
> ### Keywords: univar
> 
> ### ** Examples
> 
> data(expsy)
> cronbach(expsy[,1:10])  ## not good because item 2 is reversed (1 is high and 4 is low)
$sample.size
[1] 27

$number.of.items
[1] 10

$alpha
[1] 0.1762655

> cronbach(cbind(expsy[,c(1,3:10)],-1*expsy[,2]))  ## better
$sample.size
[1] 27

$number.of.items
[1] 10

$alpha
[1] 0.3752657

> 
> datafile <- cbind(expsy[,c(1,3:10)],-1*expsy[,2])
> library(boot)
> cronbach.boot <- function(data,x) {cronbach(data[x,])[[3]]}
> res <- boot(datafile,cronbach.boot,1000)
> quantile(res$t,c(0.025,0.975))  ## two-sided bootstrapped confidence interval of Cronbach's alpha
      2.5%      97.5% 
-0.2881438  0.6142690 
> boot.ci(res,type="bca")         ## adjusted bootstrap percentile (BCa) confidence interval (better)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = "bca")

Intervals : 
Level       BCa          
95%   (-0.1405,  0.6395 )  
Calculations and Intervals on Original Scale
> 
> 
> 
> cleanEx(); ..nameEx <- "expsy"
> 
> ### * expsy
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: expsy
> ### Title: Data set related to psychometry
> ### Aliases: expsy
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(expsy)
> expsy[1:4,]
  it1 it2 it3 it4 it5 it6 it7 it8 it9 it10 r1 rb1 r2 rb2 r3 rb3
1  NA   2   2   1   1   0   3   1   2    3 NA  NA  3   1  3   1
2   3   1   1   2   1   0   3   1   2    3  3   1  3   1  3   1
3   3   1   1   0   2   0   3   1   2    3  3   1  2   0  2   0
4   4   1   0   2   1   1   4   1   1    2  4   1  4   1  4   1
> 
> 
> 
> cleanEx(); ..nameEx <- "fpca"
> 
> ### * fpca
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fpca
> ### Title: Focused Principal Components Analysis
> ### Aliases: fpca
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(sleep)
> fpca(sleep,5,c(2:4,7:11)) 
> ## focused PCA of the duration of paradoxical sleep (dreams, 5th column)
> ## against constitutional variables in mammals (columns 2, 3, 4, 7, 8, 9, 10, 11).
> ## Variables inside the red cercle are significantly correlated
> ## to the dependent variable with p<0.05.
> ## Green variables are positively correlated to the dependent variable,
> ## yellow variables are negatively correlated.
> ## There are three clear clusters of independent variables.
> 
> corsleep <- as.data.frame(cor(sleep[,2:11],use="pairwise.complete.obs"))
> fpca(corsleep,4,c(1:3,6:10),input="Cor",sample.size=60) 
> ## when missing data are numerous, the representation of a pairwise correlation
> ## matrix may be preferred (even if mathematical properties are not so good...)
> 
> numer <- c(2:4,7:11)
> l <- length(numer)
> resu <- vector(length=l)
> for(i in 1:l)
+ {
+ int <- sleep[,numer[i]]
+ mod <- lm(sleep$Paradoxical.sleep~int)
+ resu[i] <-  summary(mod)[[4]][2,4]*sign(summary(mod)[[4]][2,1])
+ }
> fpca(sleep,5,c(2:4,7:11),pvalues=resu)
> ## A representation with p values
> ## When input="Cor" or pvalues="Yes" partial is turned to "No"
> 
> mod <- lm(sleep$Paradoxical.sleep~sleep$Body.weight+sleep$Brain.weight+
+ sleep$Slow.wave.sleep+sleep$Maximum.life.span+sleep$Gestation.time+
+ sleep$Predation+sleep$Sleep.exposure+sleep$Danger)
> resu <-  summary(mod)[[4]][2:9,4]*sign(summary(mod)[[4]][2:9,1])
> fpca(sleep,5,c(2:4,7:11),pvalues=resu)
> ## A representation with p values which come from a multiple linear model
> ## (here results are difficult to interpret)
> 
> 
> 
> cleanEx(); ..nameEx <- "icc"
> 
> ### * icc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: icc
> ### Title: Intraclass correlation coefficient (ICC)
> ### Aliases: icc
> ### Keywords: univar
> 
> ### ** Examples
> 
> data(expsy)
> icc(expsy[,c(12,14,16)])
$nb.subjects
[1] 29

$nb.raters
[1] 3

$subject.variance
[1] 0.1477833

$rater.variance
[1] 0.00410509

$residual
[1] 0.06486043

$icc.consistency
[1] 0.6949807

$icc.agreement
[1] 0.6818182

> 
> library(boot)
> icc.boot <- function(data,x) {icc(data[x,])[[7]]}
> res <- boot(expsy[,c(12,14,16)],icc.boot,1000)
> quantile(res$t,c(0.025,0.975))  # two-sided bootstrapped confidence interval of icc (agreement)
     2.5%     97.5% 
0.4193387 0.8732694 
> boot.ci(res,type="bca")         # adjusted bootstrap percentile (BCa) confidence interval (better)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = "bca")

Intervals : 
Level       BCa          
95%   ( 0.4177,  0.8722 )  
Calculations and Intervals on Original Scale
> 
> 
> 
> cleanEx(); ..nameEx <- "lkappa"
> 
> ### * lkappa
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lkappa
> ### Title: Light's kappa for n raters
> ### Aliases: lkappa
> ### Keywords: univar
> 
> ### ** Examples
> 
> data(expsy)
> lkappa(expsy[,c(11,13,15)])       # Light's kappa for non binary diagnosis
[1] 0.5471294
> lkappa(expsy[,c(12,14,16)])       # Light's kappa for binary diagnosis
[1] 0.6751938
> lkappa(expsy[,c(11,13,15)], type="weighted")      # Light's kappa for non binary ordered diagnosis
[1] 0.7571178
> 
> library(boot)
> lkappa.boot <- function(data,x) {lkappa(data[x,], type="weighted")}
> res <- boot(expsy[,c(11,13,15)],lkappa.boot,1000)
> quantile(res$t,c(0.025,0.975))    # Bootstrapped confidence interval of Light's kappa  
     2.5%     97.5% 
0.5784559 0.8547801 
> boot.ci(res,type="bca")           # adjusted bootstrap percentile (BCa) confidence interval (better)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = "bca")

Intervals : 
Level       BCa          
95%   ( 0.5853,  0.8585 )  
Calculations and Intervals on Original Scale
> 
> 
> 
> cleanEx(); ..nameEx <- "mdspca"
> 
> ### * mdspca
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mdspca
> ### Title: Graphical representation of a correlation matrix using a
> ###   Principal Component Analysis
> ### Aliases: mdspca
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(sleep)
> 
> mdspca(sleep[,c(2:5,7:11)])
> ## three consistent groups of variables, paradoxical sleep (in other words: dream)
> ## is negatively correlated with danger
> 
> mdspca(sleep[,c(2:5,7:11)],supvar=sleep[,6],namesupvar="Total.sleep",supsubj=sleep[,1],namesupsubj="",cx=0.5)
> ## Total.sleep is here a supplementary variable since it is deduced
> ## from Paradoxical.sleep and Slow.wave.sleep
> ## The variable Species is displayed in the subject plane,
> ## Rabbit and Cow have a high level of danger
> 
> 
> 
> cleanEx(); ..nameEx <- "scree.plot"
> 
> ### * scree.plot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: scree.plot
> ### Title: Screeplot of eigenvalues, simulated data are available
> ### Aliases: scree.plot
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(expsy)
> scree.plot(expsy[,1:10],simu=20,use="P") #no obvious structure with such a small sample
> 
> 
> 
> cleanEx(); ..nameEx <- "sleep"
> 
> ### * sleep
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sleep
> ### Title: Ecological and Constitutional Data in Mammals
> ### Aliases: sleep
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(sleep)
> str(sleep)
`data.frame':	62 obs. of  11 variables:
 $ Species          : Factor w/ 62 levels "African.elephant",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Body.weight      : num  6654.00    1.00    3.38    0.92 2547.00 ...
 $ Brain.weight     : num  5712.0    6.6   44.5    5.7 4603.0 ...
 $ Slow.wave.sleep  : num  NA 6.3 NA NA 2.1 9.1 15.8 5.2 10.9 8.3 ...
 $ Paradoxical.sleep: num  NA 2 NA NA 1.8 0.7 3.9 1 3.6 1.4 ...
 $ Total.sleep      : num  3.3 8.3 12.5 16.5 3.9 9.8 19.7 6.2 14.5 9.7 ...
 $ Maximum.life.span: num  38.6 4.5 14 NA 69 27 19 30.4 28 50 ...
 $ Gestation.time   : num  645 42 60 25 624 180 35 392 63 230 ...
 $ Predation        : int  3 3 1 5 3 4 1 4 1 1 ...
 $ Sleep.exposure   : int  5 1 1 2 5 4 1 5 2 1 ...
 $ Danger           : int  3 3 1 3 4 4 1 4 1 1 ...
> 
> 
> 
> cleanEx(); ..nameEx <- "sphpca"
> 
> ### * sphpca
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sphpca
> ### Title: Spherical Representation of a Correlation Matrix
> ### Aliases: sphpca
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(sleep)
> sphpca(sleep[,c(2:5,7:11)])
> ## spherical representation of ecological and constitutional correlates in mammals
> 
> sphpca(sleep[,c(2:5,7:11)],method="rscal",output=TRUE)
$stress.before.optim
[1] 0.4714109

$stress.after.optim
[1] 0.2167696

$convergence
[1] 0

$correlations
                  Body.weight Brain.weight Slow.wave.sleep Paradoxical.sleep
Body.weight             1.000        0.956          -0.394            -0.075
Brain.weight            0.956        1.000          -0.387            -0.074
Slow.wave.sleep        -0.394       -0.387           1.000             0.518
Paradoxical.sleep      -0.075       -0.074           0.518             1.000
Maximum.life.span       0.470        0.629          -0.372            -0.268
Gestation.time          0.714        0.734          -0.606            -0.409
Predation               0.096       -0.015          -0.353            -0.398
Sleep.exposure          0.406        0.323          -0.580            -0.504
Danger                  0.259        0.151          -0.535            -0.572
                  Maximum.life.span Gestation.time Predation Sleep.exposure
Body.weight                   0.470          0.714     0.096          0.406
Brain.weight                  0.629          0.734    -0.015          0.323
Slow.wave.sleep              -0.372         -0.606    -0.353         -0.580
Paradoxical.sleep            -0.268         -0.409    -0.398         -0.504
Maximum.life.span             1.000          0.646    -0.170          0.316
Gestation.time                0.646          1.000     0.091          0.573
Predation                    -0.170          0.091     1.000          0.626
Sleep.exposure                0.316          0.573     0.626          1.000
Danger                        0.015          0.306     0.927          0.790
                  Danger
Body.weight        0.259
Brain.weight       0.151
Slow.wave.sleep   -0.535
Paradoxical.sleep -0.572
Maximum.life.span  0.015
Gestation.time     0.306
Predation          0.927
Sleep.exposure     0.790
Danger             1.000

$residuals
                  Body.weight Brain.weight Slow.wave.sleep Paradoxical.sleep
Body.weight             0.000        0.058          -0.257            -0.127
Brain.weight            0.058        0.000          -0.255            -0.132
Slow.wave.sleep        -0.257       -0.255           0.000            -0.150
Paradoxical.sleep      -0.127       -0.132          -0.150             0.000
Maximum.life.span      -0.051        0.025          -0.210            -0.221
Gestation.time          0.002        0.004          -0.377            -0.306
Predation              -0.068       -0.117          -0.236            -0.255
Sleep.exposure         -0.043       -0.065          -0.350            -0.330
Danger                 -0.003       -0.051          -0.349            -0.388
                  Maximum.life.span Gestation.time Predation Sleep.exposure
Body.weight                  -0.051          0.002    -0.068         -0.043
Brain.weight                  0.025          0.004    -0.117         -0.065
Slow.wave.sleep              -0.210         -0.377    -0.236         -0.350
Paradoxical.sleep            -0.221         -0.306    -0.255         -0.330
Maximum.life.span             0.000         -0.024    -0.167          0.051
Gestation.time               -0.024          0.000    -0.120          0.032
Predation                    -0.167         -0.120     0.000          0.003
Sleep.exposure                0.051          0.032     0.003          0.000
Danger                       -0.094         -0.037     0.094          0.016
                  Danger
Body.weight       -0.003
Brain.weight      -0.051
Slow.wave.sleep   -0.349
Paradoxical.sleep -0.388
Maximum.life.span -0.094
Gestation.time    -0.037
Predation          0.094
Sleep.exposure     0.016
Danger             0.000

$mean.abs.resid
[1] 0.1407018

> ## idem, but optimizes the representation of correlations between variables with distances between points
> 
> corsleep <- as.data.frame(cor(sleep[,c(2:5,7:11)],use="pairwise.complete.obs"))
> sphpca(corsleep,input="Cor")
> sphpca(corsleep,method="rscal",input="Cor")
> ## when missing data are numerous, the representation of a pairwise correlation
> ## matrix may be preferred (even if mathematical properties are not so good...)
> 
> sphpca(corsleep,method="rscal",input="Cor",h=180,f=180,nbsphere=1,back=TRUE)
> ## other option of presentation
> 
> ##
> # library(polycor)
> # sleep$Predation <- as.ordered(sleep$Predation)
> # sleep$Sleep.exposure <- as.ordered(sleep$Sleep.exposure)
> # sleep$Danger <- as.ordered(sleep$Danger)
> # corsleeph <- as.data.frame(hetcor(sleep[,c(2:5,7:11)])$correlations)
> # sphpca(corsleeph,input="Cor",f=180)
> # sphpca(corsleeph,method="rscal",input="Cor",f=180)
> ## --> Correlations between discrete variables may appear shoking to some statisticians (?)
> ## --> Representation of polychoric/polyserial correlations could be prefered in this situation
> 
> 
> 
> cleanEx(); ..nameEx <- "wkappa"
> 
> ### * wkappa
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: wkappa
> ### Title: weighted Kappa for 2 raters
> ### Aliases: wkappa
> ### Keywords: univar
> 
> ### ** Examples
> 
> data(expsy)
> wkappa(expsy[,c(11,13)])        # weighted kappa (squared weights)
$table
  1 2 3 4
1 1 0 0 0
2 1 4 0 0
3 0 4 9 2
4 0 0 2 6

$weights
[1] "squared"

$kappa
[1] 0.7819549

> 
> library(boot)
> wkappa.boot <- function(data,x) {wkappa(data[x,])[[3]]}
> res <- boot(expsy[,c(11,13)],wkappa.boot,1000)
> quantile(res$t,c(0.025,0.975))  # two-sided bootstrapped confidence interval of weighted kappa
     2.5%     97.5% 
0.5805037 0.8965670 
> boot.ci(res,type="bca")         # adjusted bootstrap percentile (BCa) confidence interval (better)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = "bca")

Intervals : 
Level       BCa          
95%   ( 0.6055,  0.9109 )  
Calculations and Intervals on Original Scale
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  73.72 2.29 76.26 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')