R version 2.10.0 Under development (unstable) (2009-05-22 r48594)
Copyright (C) 2009 The R Foundation for Statistical Computing
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.

  Natural language support but running in an English locale

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.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign("nameEx",
+        local({
+ 	   s <- "__{must remake R-ex/*.R}__"
+            function(new) {
+                if(!missing(new)) s <<- new else s
+            }
+        }),
+        pos = "CheckExEnv")
> ## Add some hooks to label plot pages for base and grid graphics
> assign("base_plot_hook",
+        function() {
+            pp <- par(c("mfg","mfcol","oma","mar"))
+            if(all(pp$mfg[1:2] == c(1, pp$mfcol[2]))) {
+                outer <- (oma4 <- pp$oma[4]) > 0; mar4 <- pp$mar[4]
+                mtext(sprintf("help(\"%s\")", nameEx()), side = 4,
+                      line = if(outer)max(1, oma4 - 1) else min(1, mar4 - 1),
+                outer = outer, adj = 1, cex = .8, col = "orchid", las=3)
+            }
+        },
+        pos = "CheckExEnv")
> assign("grid_plot_hook",
+        function() {
+            grid::pushViewport(grid::viewport(width=grid::unit(1, "npc") -
+                               grid::unit(1, "lines"), x=0, just="left"))
+            grid::grid.text(sprintf("help(\"%s\")", nameEx()),
+                            x=grid::unit(1, "npc") + grid::unit(0.5, "lines"),
+                            y=grid::unit(0.8, "npc"), rot=90,
+                            gp=grid::gpar(col="orchid"))
+        },
+        pos = "CheckExEnv")
> setHook("plot.new",     get("base_plot_hook", pos = "CheckExEnv"))
> setHook("persp",        get("base_plot_hook", pos = "CheckExEnv"))
> setHook("grid.newpage", get("grid_plot_hook", pos = "CheckExEnv"))
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   .CheckExEnv <- as.environment("CheckExEnv")
+ 	   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")
+        },
+        pos = "CheckExEnv")
> assign("ptime", proc.time(), pos = "CheckExEnv")
> ## at least one package changes these via ps.options(), so do this
> ## before loading the package.
> ## Use postscript as incomplete files may be viewable, unlike PDF.
> ## Choose a size that is close to on-screen devices, fix paper
> grDevices::ps.options(width = 7, height = 7, paper = "a4", reset = TRUE)
> grDevices::postscript("psych-Ex.ps")
> 
> assign("par.postscript", graphics::par(no.readonly = TRUE), pos = "CheckExEnv")
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)
> library('psych')
> 
> assign(".oldSearch", search(), pos = 'CheckExEnv')
> assign(".oldNS", loadedNamespaces(), pos = 'CheckExEnv')
> cleanEx(); nameEx("00.psych-package")
> ### * 00.psych-package
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: 00.psych
> ### Title: A package for personality, psychometric, and psychological
> ###   research
> ### Aliases: psych psych-package 00.psych-package
> ### Keywords: package multivariate models cluster
> 
> ### ** Examples
> 
> #See the separate man pages 
> test.psych()
too many factors requested for this number of variables to use SMC, 1s used instead
Loading required package: GPArotation
In omega, 3 factors are too many for 4 variables using mle.  Using minres instead
In fa, too many factors requested for this number of variables to use SMC for communality estimates, 1s are used instead
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
In fa, too many factors requested for this number of variables to use SMC for communality estimates, 1s are used instead
In smc, the correlation matrix was not invertible, smc's returned as 1s
Warning in log(det(m.inv.r)) : NaNs produced
In factor.stats, the correlation matrix is singular, an approximation is used
In factor.scores, the correlation matrix is singular, an approximation is used
Warning in log(det(m.inv.r)) : NaNs produced
In factor.stats, the correlation matrix is singular, an approximation is used
In factor.scores, the correlation matrix is singular, an approximation is used
Loading required package: MASS
> 
> 
> 
> cleanEx(); nameEx("Harman")
> ### * Harman
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Harman
> ### Title: Two data sets from Harman (1967). 9 cognitive variables from
> ###   Holzinger and 8 emotional variables from Burt
> ### Aliases: Harman Harman.Burt Harman.Holzinger Burt
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(Harman)
> cor.plot(Harman.Holzinger)
> cor.plot(Harman.Burt)  
> smc(Harman.Burt)  #note how this produces impossible results
Sociability      Sorrow  Tenderness         Joy      Wonder     Disgust 
  1.0606356   1.1388831   1.0464691   0.2487892   2.4397325   1.8048109 
      Anger        Fear 
  1.9187862   0.3152856 
> 
> 
> 
> cleanEx(); nameEx("ICC")
> ### * ICC
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ICC
> ### Title: Intraclass Correlations (ICC1, ICC2, ICC3 from Shrout and
> ###   Fleiss)
> ### Aliases: ICC
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> sf <- matrix(c(9,    2,   5,    8,
+ 6,    1,   3,    2,
+ 8,    4,   6,    8,
+ 7,    1,   2,    6,
+ 10,   5,   6,    9,
+ 6,   2,   4,    7),ncol=4,byrow=TRUE)
> colnames(sf) <- paste("J",1:4,sep="")
> rownames(sf) <- paste("S",1:6,sep="")
> sf  #example from Shrout and Fleiss (1979)
   J1 J2 J3 J4
S1  9  2  5  8
S2  6  1  3  2
S3  8  4  6  8
S4  7  1  2  6
S5 10  5  6  9
S6  6  2  4  7
> ICC(sf)
                         type  ICC     F df1 df2    p lower bound upper bound
Single_raters_absolute   ICC1 0.17  1.79   5  18 0.16       -0.13        0.72
Single_random_raters     ICC2 0.29 11.03   5  15 0.00        0.02        0.76
Single_fixed_raters      ICC3 0.71 11.03   5  15 0.00        0.34        0.95
Average_raters_absolute ICC1k 0.44  1.79   5  18 0.16       -0.88        0.91
Average_random_raters   ICC2k 0.62 11.03   5  15 0.00        0.07        0.93
Average_fixed_raters    ICC3k 0.91 11.03   5  15 0.00        0.68        0.99
> 
> 
> 
> cleanEx(); nameEx("ICLUST")
> ### * ICLUST
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ICLUST
> ### Title: ICLUST: Item Cluster Analysis - Hierarchical cluster analysis
> ###   using psychometric principles
> ### Aliases: ICLUST iclust
> ### Keywords: multivariate cluster
> 
> ### ** Examples
> 
> test.data <- Harman74.cor$cov
> ic.out <- ICLUST(test.data)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
> summary(ic.out)
ICLUST (Item Cluster Analysis)Call: ICLUST(r.mat = test.data)
ICLUST 

Purified Alpha:
[1] 0.91

 Guttman Lambda6  * 
[1] 0.94

Original Beta:
[1] 0.78

Cluster size:
[1] 24

Purified scale intercorrelations
 reliabilities on diagonal
 correlations corrected for attenuation above diagonal: 
     [,1]
[1,] 0.91
> ic.out <- ICLUST(test.data,nclusters =4)  #use all defaults and stop at 4 clusters
> ic.out1 <- ICLUST(test.data,beta=3,beta.size=3)  #use more stringent criteria
> print(ic.out1)
ICLUST (Item Cluster Analysis)
Call: ICLUST(r.mat = test.data, beta = 3, beta.size = 3)

Purified Alpha:
 C18  C13  C15  C20 
0.82 0.90 0.81 0.74 

G6* reliability:
 C18  C13  C20  C15 
0.84 0.90 0.76 0.84 

Original Beta:
 C18  C13  C20  C15 
0.74 0.79 0.66 0.74 

Cluster size:
C18 C13 C15 C20 
  8   5   5   6 

Item by Cluster Structure matrix:
                       C18  C13  C20  C15 
VisualPerception       0.68 0.40 0.46 0.43
Cubes                  0.48 0.26 0.25 0.23
PaperFormBoard         0.53 0.33 0.35 0.20
Flags                  0.58 0.40 0.32 0.27
GeneralInformation     0.55 0.81 0.39 0.47
PargraphComprehension  0.55 0.81 0.44 0.40
SentenceCompletion     0.53 0.84 0.35 0.42
WordClassification     0.60 0.70 0.44 0.47
WordMeaning            0.56 0.84 0.44 0.36
Addition               0.25 0.30 0.36 0.75
Code                   0.37 0.38 0.55 0.69
CountingDots           0.39 0.22 0.34 0.71
StraightCurvedCapitals 0.56 0.42 0.40 0.68
WordRecognition        0.30 0.33 0.55 0.33
NumberRecognition      0.32 0.27 0.54 0.25
FigureRecognition      0.54 0.30 0.61 0.32
ObjectNumber           0.32 0.31 0.64 0.42
NumberFigure           0.48 0.26 0.62 0.50
FigureWord             0.41 0.30 0.49 0.35
Deduction              0.65 0.54 0.49 0.35
NumericalPuzzles       0.60 0.38 0.47 0.59
ProblemReasoning       0.63 0.51 0.49 0.38
SeriesCompletion       0.73 0.57 0.50 0.48
ArithmeticProblems     0.53 0.52 0.53 0.65

With eigenvalues of:
C18 C13 C20 C15 
6.5 6.0 5.3 5.4 

Purified scale intercorrelations
 reliabilities on diagonal
 correlations corrected for attenuation above diagonal: 
     C18  C13  C15  C20
C18 0.82 0.71 0.62 0.71
C13 0.61 0.90 0.54 0.52
C15 0.51 0.46 0.81 0.65
C20 0.55 0.43 0.51 0.74
> plot(ic.out)    #this shows the spatial representation
Use ICLUST.graph to see the  hierarchical structure
> ic.no.graph <- ICLUST(test.data,plot=FALSE)
> dot.graph <- ICLUST.graph(ic.no.graph,out.file="test.ICLUST.graph.dot")   #use a dot graphics viewer
> 
> 
> 
> 
> cleanEx(); nameEx("ICLUST.graph")
> ### * ICLUST.graph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ICLUST.graph
> ### Title: create control code for ICLUST graphical output
> ### Aliases: ICLUST.graph
> ### Keywords: multivariate cluster hplot
> 
> ### ** Examples
> 
> ## Not run: 
> ##D test.data <- Harman74.cor$cov
> ##D ic.out <- ICLUST(test.data)
> ##D out.file <- file.choose(new=TRUE)   #create a new file to write the plot commands to 
> ##D ICLUST.graph(ic.out,out.file)   
> ##D now go to graphviz (outside of R) and open the out.file you created
> ##D print(ic.out,digits=2)
> ## End(Not run)
> 
>  
> #test.data <- Harman74.cor$cov 
> #my.iclust <- ICLUST(test.data)
> #ICLUST.graph(my.iclust)
> #
> #
> #digraph ICLUST {
> #  rankdir=RL;
> #  size="8,8";
> #  node [fontname="Helvetica" fontsize=14 shape=box, width=2];
> #  edge [fontname="Helvetica" fontsize=12];
> # label = "ICLUST";
> #       fontsize=20;
> #V1  [label = VisualPerception];
> #V2  [label = Cubes];
> #V3  [label = PaperFormBoard];
> #V4  [label = Flags];
> #V5  [label = GeneralInformation];
> #V6  [label = PargraphComprehension];
> #V7  [label = SentenceCompletion];
> #V8  [label = WordClassification];
> #V9  [label = WordMeaning];
> #V10  [label = Addition];
> #V11  [label = Code];
> #V12  [label = CountingDots];
> #V13  [label = StraightCurvedCapitals];
> #V14  [label = WordRecognition];
> #V15  [label = NumberRecognition];
> #V16  [label = FigureRecognition];
> #V17  [label = ObjectNumber];
> #V18  [label = NumberFigure];
> #V19  [label = FigureWord];
> #V20  [label = Deduction];
> #V21  [label = NumericalPuzzles];
> #V22  [label = ProblemReasoning];
> #V23  [label = SeriesCompletion];
> #V24  [label = ArithmeticProblems];
> #node [shape=ellipse, width ="1"];
> #C1-> V9 [ label = 0.78 ];
> #C1-> V5 [ label = 0.78 ];
> #C2-> V12 [ label = 0.66 ];
> #C2-> V10 [ label = 0.66 ];
> #C3-> V18 [ label = 0.53 ];
> #C3-> V17 [ label = 0.53 ];
> #C4-> V23 [ label = 0.59 ];
> #C4-> V20 [ label = 0.59 ];
> #C5-> V13 [ label = 0.61 ];
> #C5-> V11 [ label = 0.61 ];
> #C6-> V7 [ label = 0.78 ];
> #C6-> V6 [ label = 0.78 ];
> #C7-> V4 [ label = 0.55 ];
> #C7-> V1 [ label = 0.55 ];
> #C8-> V16 [ label = 0.5 ];
> #C8-> V14 [ label = 0.49 ];
> #C9-> C1 [ label = 0.86 ];
> #C9-> C6 [ label = 0.86 ];
> #C10-> C4 [ label = 0.71 ];
> #C10-> V22 [ label = 0.62 ];
> #C11-> V21 [ label = 0.56 ];
> #C11-> V24 [ label = 0.58 ];
> #C12-> C10 [ label = 0.76 ];
> #C12-> C11 [ label = 0.67 ];
> #C13-> C8 [ label = 0.61 ];
> #C13-> V15 [ label = 0.49 ];
> #C14-> C2 [ label = 0.74 ];
> #C14-> C5 [ label = 0.72 ];
> #C15-> V3 [ label = 0.48 ];
> #C15-> C7 [ label = 0.65 ];
> #C16-> V19 [ label = 0.48 ];
> #C16-> C3 [ label = 0.64 ];
> #C17-> V8 [ label = 0.62 ];
> #C17-> C12 [ label = 0.8 ];
> #C18-> C17 [ label = 0.82 ];
> #C18-> C15 [ label = 0.68 ];
> #C19-> C16 [ label = 0.66 ];
> #C19-> C13 [ label = 0.65 ];
> #C20-> C19 [ label = 0.72 ];
> #C20-> C18 [ label = 0.83 ];
> #C21-> C20 [ label = 0.87 ];
> #C21-> C9 [ label = 0.76 ];
> #C22-> 0 [ label = 0 ];
> #C22-> 0 [ label = 0 ];
> #C23-> 0 [ label = 0 ];
> #C23-> 0 [ label = 0 ];
> #C1  [label =   "C1\n  alpha= 0.84\n beta=  0.84\nN= 2"] ;
> #C2  [label =   "C2\n  alpha= 0.74\n beta=  0.74\nN= 2"] ;
> #C3  [label =   "C3\n  alpha= 0.62\n beta=  0.62\nN= 2"] ;
> #C4  [label =   "C4\n  alpha= 0.67\n beta=  0.67\nN= 2"] ;
> #C5  [label =   "C5\n  alpha= 0.7\n beta=  0.7\nN= 2"] ;
> #C6  [label =   "C6\n  alpha= 0.84\n beta=  0.84\nN= 2"] ;
> #C7  [label =   "C7\n  alpha= 0.64\n beta=  0.64\nN= 2"] ;
> #C8  [label =   "C8\n  alpha= 0.58\n beta=  0.58\nN= 2"] ;
> #C9  [label =   "C9\n  alpha= 0.9\n beta=  0.87\nN= 4"] ;
> #C10  [label =   "C10\n  alpha= 0.74\n beta=  0.71\nN= 3"] ;
> #C11  [label =   "C11\n  alpha= 0.62\n beta=  0.62\nN= 2"] ;
> #C12  [label =   "C12\n  alpha= 0.79\n beta=  0.74\nN= 5"] ;
> #C13  [label =   "C13\n  alpha= 0.64\n beta=  0.59\nN= 3"] ;
> #C14  [label =   "C14\n  alpha= 0.79\n beta=  0.74\nN= 4"] ;
> #C15  [label =   "C15\n  alpha= 0.66\n beta=  0.58\nN= 3"] ;
> #C16  [label =   "C16\n  alpha= 0.65\n beta=  0.57\nN= 3"] ;
> #C17  [label =   "C17\n  alpha= 0.81\n beta=  0.71\nN= 6"] ;
> #C18  [label =   "C18\n  alpha= 0.84\n beta=  0.75\nN= 9"] ;
> #C19  [label =   "C19\n  alpha= 0.74\n beta=  0.65\nN= 6"] ;
> #C20  [label =   "C20\n  alpha= 0.87\n beta=  0.74\nN= 15"] ;
> #C21  [label =   "C21\n  alpha= 0.9\n beta=  0.77\nN= 19"] ;
> #C22  [label =   "C22\n  alpha= 0\n beta=  0\nN= 0"] ;
> #C23  [label =   "C23\n  alpha= 0\n beta=  0\nN= 0"] ;
> #{ rank=same;
> #V1;V2;V3;V4;V5;V6;V7;V8;V9;V10;V11;V12;V13;V14;V15;V16;V17;V18;V19;V20;V21;V22;V23;V24;}}
> #
> #copy the above output to Graphviz and draw it
> #see \url{http://personality-project.org/r/r.ICLUST.html} for an example.
> 
> 
> 
> 
> cleanEx(); nameEx("ICLUST.rgraph")
> ### * ICLUST.rgraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ICLUST.rgraph
> ### Title: Draw an ICLUST graph using the Rgraphviz package
> ### Aliases: ICLUST.rgraph
> ### Keywords: multivariate cluster hplot
> 
> ### ** Examples
> 
> test.data <- Harman74.cor$cov
> 
> if(require(Rgraphviz) ) {ic.out <- ICLUST(test.data) }  
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
> 
> 
> 
> 
> cleanEx(); nameEx("Promax")
> ### * Promax
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Promax
> ### Title: Perform promax or targeted rotations and return the inter factor
> ###   angles
> ### Aliases: Promax target.rot
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> jen <- sim.hierarchical()
Loading required package: MASS
> f3 <- factor.minres(jen,3)
> Promax(f3)

Call: NULL
   V  MR1  MR2  MR3   h2   u2
V1 1 0.76           0.58 0.42
V2 2 0.66           0.44 0.56
V3 3 0.57           0.32 0.68
V4 4      0.69      0.47 0.53
V5 5      0.59      0.35 0.65
V6 6      0.49      0.24 0.76
V7 7           0.60 0.36 0.64
V8 8           0.50 0.25 0.75
V9 9           0.40 0.16 0.84

                MR1  MR2  MR3
SS loadings    1.34 1.07 0.77
Proportion Var 0.15 0.12 0.09
Cumulative Var 0.15 0.27 0.35

 With factor correlations of 
     MR1  MR2  MR3
MR1 1.00 0.68 0.60
MR2 0.68 1.00 0.55
MR3 0.60 0.55 1.00
> target.rot(f3)

Call: NULL
   V MR1 MR2 MR3   h2   u2
V1 1 0.8         0.64 0.36
V2 2 0.7         0.49 0.51
V3 3 0.6         0.36 0.64
V4 4     0.7     0.49 0.51
V5 5     0.6     0.36 0.64
V6 6     0.5     0.25 0.75
V7 7         0.6 0.36 0.64
V8 8         0.5 0.25 0.75
V9 9         0.4 0.16 0.84

                MR1  MR2  MR3
SS loadings    1.49 1.10 0.77
Proportion Var 0.17 0.12 0.09
Cumulative Var 0.17 0.29 0.37

 With factor correlations of 
     MR1  MR2  MR3
MR1 1.00 0.72 0.63
MR2 0.72 1.00 0.56
MR3 0.63 0.56 1.00
> m3 <- factanal(covmat=jen,factors=3)
> Promax(m3)  #example of taking the output from factanal

Call: NULL
   V Factor1 Factor2 Factor3   h2   u2
V1 1    0.76                 0.58 0.42
V2 2    0.66                 0.44 0.56
V3 3    0.57                 0.32 0.68
V4 4            0.69         0.47 0.53
V5 5            0.59         0.35 0.65
V6 6            0.49         0.24 0.76
V7 7                    0.60 0.36 0.64
V8 8                    0.50 0.25 0.75
V9 9                    0.40 0.16 0.84

               Factor1 Factor2 Factor3
SS loadings       1.34    1.07    0.77
Proportion Var    0.15    0.12    0.09
Cumulative Var    0.15    0.27    0.35

 With factor correlations of 
        Factor1 Factor2 Factor3
Factor1    1.00    0.68    0.60
Factor2    0.68    1.00    0.55
Factor3    0.60    0.55    1.00
> #compare this rotation with the solution from a targeted rotation aimed for an independent cluster solution
> target.rot(m3)

Call: NULL
   V Factor1 Factor2 Factor3   h2   u2
V1 1     0.8                 0.64 0.36
V2 2     0.7                 0.49 0.51
V3 3     0.6                 0.36 0.64
V4 4             0.7         0.49 0.51
V5 5             0.6         0.36 0.64
V6 6             0.5         0.25 0.75
V7 7                     0.6 0.36 0.64
V8 8                     0.5 0.25 0.75
V9 9                     0.4 0.16 0.84

               Factor1 Factor2 Factor3
SS loadings       1.49    1.10    0.77
Proportion Var    0.17    0.12    0.09
Cumulative Var    0.17    0.29    0.37

 With factor correlations of 
        Factor1 Factor2 Factor3
Factor1    1.00    0.72    0.63
Factor2    0.72    1.00    0.56
Factor3    0.63    0.56    1.00
> 
> 
> 
> cleanEx(); nameEx("SD")
> ### * SD
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: SD
> ### Title: Find the Standard deviation for a vector, matrix, or data.frame
> ###   - do not return error if there are no cases
> ### Aliases: SD
> ### Keywords: models
> 
> ### ** Examples
> 
> data(attitude)
> sd(attitude) #all complete
    rating complaints privileges   learning     raises   critical    advance 
 12.172562  13.314757  12.235430  11.737013  10.397226   9.894908  10.288706 
> attitude[,1] <- NA
> SD(attitude) #missing a column
    rating complaints privileges   learning     raises   critical    advance 
        NA  13.314757  12.235430  11.737013  10.397226   9.894908  10.288706 
> describe(attitude)
Warning in min(x[, i], na.rm = na.rm) :
  no non-missing arguments to min; returning Inf
Warning in max(x[, i], na.rm = na.rm) :
  no non-missing arguments to max; returning -Inf
           var  n  mean    sd median trimmed   mad min  max range  skew
rating*      1  0   NaN    NA     NA     NaN    NA Inf -Inf  -Inf    NA
complaints   2 30 66.60 13.31   65.0   67.08 14.83  37   90    53 -0.22
privileges   3 30 53.13 12.24   51.5   52.75 10.38  30   83    53  0.38
learning     4 30 56.37 11.74   56.5   56.58 14.83  34   75    41 -0.05
raises       5 30 64.63 10.40   63.5   64.50 11.12  43   88    45  0.20
critical     6 30 74.77  9.89   77.5   75.83  7.41  49   92    43 -0.87
advance      7 30 42.93 10.29   41.0   41.83  8.90  25   72    47  0.85
           kurtosis   se
rating*          NA   NA
complaints    -0.68 2.43
privileges    -0.41 2.23
learning      -1.22 2.14
raises        -0.60 1.90
critical       0.17 1.81
advance        0.47 1.88
> 
> 
> 
> cleanEx(); nameEx("VSS")
> ### * VSS
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: VSS
> ### Title: Apply the Very Simple Structure and MAP criteria to determine
> ###   the appropriate number of factors.
> ### Aliases: VSS MAP
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> test.data <- Harman74.cor$cov
> my.vss <- VSS(test.data,title="VSS of 24 mental tests")      
n.obs was not specified and was arbitrarily set to 1000.  This only affects the chi square values.
> #print(my.vss[,1:12],digits =2) 
> #VSS.plot(my.vss, title="VSS of 24 mental tests")
> 
> #now, some simulated data with two factors
> VSS(sim.circ(nvar=24),fm="mle" ,title="VSS of 24 circumplex variables")

Very Simple Structure of  VSS of 24 circumplex variables 
Call: VSS(x = sim.circ(nvar = 24), fm = "mle", title = "VSS of 24 circumplex variables")
VSS complexity 1 achieves a maximimum of 0.61  with  3  factors
VSS complexity 2 achieves a maximimum of 0.85  with  3  factors

The Velicer MAP criterion achieves a minimum of 0.05  with  2  factors
 
Velicer MAP
[1] 0.05 0.01 0.01 0.01 0.01 0.02 0.02 0.02

Very Simple Structure Complexity 1
[1] 0.46 0.61 0.61 0.58 0.56 0.55 0.51 0.48

Very Simple Structure Complexity 2
[1] 0.00 0.85 0.85 0.83 0.80 0.78 0.76 0.74
> VSS(sim.item(nvar=24),fm="mle" ,title="VSS of 24 simple structure variables")

Very Simple Structure of  VSS of 24 simple structure variables 
Call: VSS(x = sim.item(nvar = 24), fm = "mle", title = "VSS of 24 simple structure variables")
VSS complexity 1 achieves a maximimum of 0.79  with  2  factors
VSS complexity 2 achieves a maximimum of 0.83  with  8  factors

The Velicer MAP criterion achieves a minimum of 0.04  with  2  factors
 
Velicer MAP
[1] 0.04 0.01 0.01 0.01 0.01 0.02 0.02 0.02

Very Simple Structure Complexity 1
[1] 0.43 0.79 0.79 0.79 0.74 0.70 0.67 0.64

Very Simple Structure Complexity 2
[1] 0.00 0.80 0.80 0.81 0.81 0.82 0.82 0.83
> 
> 
> 
> cleanEx(); nameEx("VSS.parallel")
> ### * VSS.parallel
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: VSS.parallel
> ### Title: Compare real and random VSS solutions
> ### Aliases: VSS.parallel
> ### Keywords: models models
> 
> ### ** Examples
> 
> #VSS.plot(VSS.parallel(200,24))
> 
> 
> 
> cleanEx(); nameEx("VSS.plot")
> ### * VSS.plot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: VSS.plot
> ### Title: Plot VSS fits
> ### Aliases: VSS.plot
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> test.data <- Harman74.cor$cov
> my.vss <- VSS(test.data)         #suggests that 4 factor complexity two solution is optimal
n.obs was not specified and was arbitrarily set to 1000.  This only affects the chi square values.
> VSS.plot(my.vss,title="VSS of Holzinger-Harmon problem")                 #see the graphics window
> 
> 
> 
> 
> cleanEx(); nameEx("VSS.scree")
> ### * VSS.scree
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: VSS.scree
> ### Title: Plot a scree test
> ### Aliases: VSS.scree
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> #VSS.scree(attitude)
> #VSS.scree(cor(attitude)
> 
> 
> 
> 
> cleanEx(); nameEx("Yule")
> ### * Yule
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Yule
> ### Title: From a two by two table, find the Yule coefficients of
> ###   association, convert to phi, or polychoric, recreate table the table
> ###   to create the Yule coefficient.
> ### Aliases: Yule Yule.inv Yule2phi Yule2poly
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> Nach <- matrix(c(40,10,20,50),ncol=2,byrow=TRUE)
> Yule(Nach)
[1] 0.8181818
> Yule.inv(.81818,c(50,70,60,60))
     [,1] [,2]
[1,]   40   10
[2,]   20   50
> Yule2phi(.81818,c(50,70,60,60))
[1] 0.5070926
> if(require(polycor)) Yule2poly(.81818,c(50,70,60,60))
Loading required package: polycor
Loading required package: mvtnorm
Loading required package: sfsmisc
Warning in fun(...) :
  no valid postscript previewer found; consider setting
  options("eps_view"=  "....")	    yourself
[1] 0.7230661
> phi(Nach)  #much less
[1] 0.51
> 
> 
> 
> cleanEx(); nameEx("alpha")
> ### * alpha
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: alpha
> ### Title: Find two estimates of reliability: Cronbach's alpha and
> ###   Guttman's Lambda 6.
> ### Aliases: alpha alpha.scale
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> r4 <- sim.congeneric()
> alpha(r4)

Reliability analysis   
Call: alpha(x = r4)

  raw_alpha std.alpha G6(smc) average_r
      0.74      0.74     0.7      0.42

 Reliability if an item is dropped:
   raw_alpha std.alpha G6(smc) average_r
V1      0.62      0.62    0.53      0.36
V2      0.66      0.66    0.57      0.39
V3      0.70      0.70    0.62      0.44
V4      0.74      0.74    0.66      0.49

 Item statistics 
      r r.cor
V1 0.81  0.74
V2 0.78  0.67
V3 0.73  0.59
V4 0.68  0.50
> r9 <- sim.hierarchical()
Loading required package: MASS
> alpha(r9)

Reliability analysis   
Call: alpha(x = r9)

  raw_alpha std.alpha G6(smc) average_r
      0.76      0.76    0.76      0.26

 Reliability if an item is dropped:
   raw_alpha std.alpha G6(smc) average_r
V1      0.71      0.71    0.70      0.24
V2      0.72      0.72    0.71      0.25
V3      0.74      0.74    0.73      0.26
V4      0.73      0.73    0.72      0.25
V5      0.74      0.74    0.73      0.26
V6      0.75      0.75    0.74      0.27
V7      0.75      0.75    0.74      0.27
V8      0.76      0.76    0.75      0.28
V9      0.77      0.77    0.76      0.29

 Item statistics 
      r r.cor
V1 0.72  0.71
V2 0.67  0.63
V3 0.61  0.55
V4 0.65  0.59
V5 0.59  0.52
V6 0.53  0.43
V7 0.56  0.46
V8 0.50  0.39
V9 0.45  0.32
> #an example of two independent factors that produce reasonable alphas
> #this is a case where alpha is a poor indicator of unidimensionality
> two.f <- sim.item(8)
> alpha(two.f,keys=c(rep(1,4),rep(-1,4))) 

Reliability analysis   
Call: alpha(x = two.f, keys = c(rep(1, 4), rep(-1, 4)))

  raw_alpha std.alpha G6(smc) average_r   mean  sd
       0.6       0.6    0.65      0.16 -0.015 2.3

 Reliability if an item is dropped:
  raw_alpha std.alpha G6(smc) average_r
1      0.56      0.56    0.62      0.16
2      0.56      0.57    0.61      0.16
3      0.57      0.57    0.62      0.16
4      0.57      0.57    0.61      0.16
5      0.57      0.57    0.62      0.16
6      0.58      0.58    0.62      0.17
7      0.57      0.57    0.61      0.16
8      0.57      0.57    0.61      0.16

 Item statistics 
    n    r r.cor    mean   sd
1 500 0.53  0.43 -0.0512 1.08
2 500 0.53  0.43 -0.0203 1.04
3 500 0.51  0.42  0.0154 1.05
4 500 0.52  0.43  0.0701 1.06
5 500 0.50  0.39 -0.0176 1.01
6 500 0.47  0.37  0.0225 1.02
7 500 0.52  0.42 -0.0252 1.00
8 500 0.52  0.43 -0.0086 0.96
> 
> 
> 
> 
> cleanEx(); nameEx("bfi")
> ### * bfi
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bfi
> ### Title: 25 Personality items representing 5 factors
> ### Aliases: bfi
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(bfi)
> describe(bfi)
   var    n mean   sd median trimmed  mad min max range  skew kurtosis   se
A1   1 1000 2.29 1.28      2    2.13 1.48   1   6     5  0.87    -0.11 0.04
A2   2  994 4.82 1.12      5    4.98 1.48   1   6     5 -1.08     1.03 0.04
A3   3  989 4.60 1.21      5    4.76 1.48   1   6     5 -1.00     0.60 0.04
A4   4  993 4.76 1.40      5    5.00 1.48   1   6     5 -1.08     0.29 0.04
A5   5  988 4.58 1.15      5    4.70 1.48   1   6     5 -0.76     0.24 0.04
C1   6  997 4.39 1.22      5    4.50 1.48   1   6     5 -0.72     0.00 0.04
C2   7  997 4.22 1.29      4    4.32 1.48   1   6     5 -0.64    -0.18 0.04
C3   8  995 4.28 1.26      5    4.38 1.48   1   6     5 -0.68    -0.13 0.04
C4   9  986 2.63 1.36      2    2.51 1.48   1   6     5  0.51    -0.76 0.04
C5  10  997 3.47 1.52      4    3.48 1.48   1   6     5 -0.09    -1.06 0.05
E1  11  996 2.93 1.58      3    2.84 1.48   1   6     5  0.36    -1.07 0.05
E2  12  994 3.32 1.58      3    3.29 1.48   1   6     5  0.05    -1.16 0.05
E3  13  994 3.98 1.29      4    4.02 1.48   1   6     5 -0.35    -0.55 0.04
E4  14  997 4.42 1.43      5    4.59 1.48   1   6     5 -0.86    -0.19 0.05
E5  15  991 4.38 1.24      5    4.50 1.48   1   6     5 -0.72    -0.03 0.04
N1  16  990 2.83 1.51      3    2.72 1.48   1   6     5  0.46    -0.88 0.05
N2  17  990 3.51 1.53      4    3.52 1.48   1   6     5  0.01    -1.07 0.05
N3  18  997 3.20 1.52      3    3.15 1.48   1   6     5  0.22    -1.06 0.05
N4  19  996 3.06 1.54      3    2.97 1.48   1   6     5  0.35    -0.96 0.05
N5  20  992 2.94 1.55      3    2.84 1.48   1   6     5  0.42    -0.95 0.05
O1  21  994 4.74 1.18      5    4.88 1.48   1   6     5 -0.82     0.13 0.04
O2  22  994 3.93 1.33      4    3.95 1.48   1   6     5 -0.15    -0.75 0.04
O3  23  991 4.32 1.22      5    4.42 1.48   1   6     5 -0.66    -0.05 0.04
O4  24  992 4.83 1.18      5    5.00 1.48   1   6     5 -1.01     0.58 0.04
O5  25  996 2.51 1.22      2    2.39 1.48   1   6     5  0.65    -0.22 0.04
>  data(bfi)
>  keys.list <- list(Agree=c(-1,2:5),Conscientious=c(6:8,-9,-10),Extraversion=c(-11,-12,13:15),Neuroticism=c(16:20),Openness = c(21,-22,23,24,-25))
>  keys <- make.keys(25,keys.list,item.labels=colnames(bfi))
>  scores <- score.items(keys,bfi,short=TRUE)
>  scores
Call: score.items(keys = keys, items = bfi, short = TRUE)

(Unstandardized) Alpha:
      Agree Conscientious Extraversion Neuroticism Openness
alpha  0.65          0.74         0.77         0.8     0.47

Average item correlation:
          Agree Conscientious Extraversion Neuroticism Openness
average.r  0.27          0.36         0.41        0.46     0.15

 Guttman 6* reliability: 
         Agree Conscientious Extraversion Neuroticism Openness
Lambda.6  0.66          0.74         0.78        0.82     0.53

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
              Agree Conscientious Extraversion Neuroticism Openness
Agree          0.65          0.36         0.62       -0.30     0.33
Conscientious  0.25          0.74         0.28       -0.22     0.24
Extraversion   0.44          0.21         0.77       -0.24     0.57
Neuroticism   -0.22         -0.17        -0.19        0.81     0.11
Openness       0.18          0.14         0.34        0.07     0.47

Item by scale correlations:
 corrected for item overlap and scale reliability
   Agree Conscientious Extraversion Neuroticism Openness
A1 -0.40         -0.12        -0.07        0.16    -0.16
A2  0.63          0.20         0.41       -0.09     0.26
A3  0.63          0.20         0.41       -0.13     0.30
A4  0.39          0.21         0.22       -0.17    -0.05
A5  0.61          0.19         0.54       -0.23     0.17
C1  0.12          0.59         0.11       -0.03     0.23
C2  0.19          0.56         0.08        0.07     0.14
C3  0.24          0.61         0.23       -0.03     0.15
C4 -0.29         -0.70        -0.22        0.30    -0.17
C5 -0.22         -0.57        -0.20        0.32    -0.01
E1 -0.36         -0.07        -0.65        0.09    -0.42
E2 -0.35         -0.21        -0.69        0.28    -0.23
E3  0.44          0.19         0.61       -0.09     0.50
E4  0.53          0.12         0.68       -0.23     0.10
E5  0.27          0.34         0.58       -0.05     0.49
N1 -0.26         -0.17        -0.06        0.77     0.13
N2 -0.26         -0.11        -0.06        0.74     0.09
N3 -0.19         -0.16        -0.10        0.71     0.10
N4 -0.24         -0.27        -0.37        0.63     0.08
N5 -0.06         -0.04        -0.23        0.54    -0.03
O1  0.21          0.17         0.30       -0.01     0.54
O2  0.13          0.13        -0.22       -0.06    -0.09
O3  0.29          0.21         0.43        0.01     0.63
O4  0.17          0.10         0.00        0.22     0.42
O5 -0.12         -0.14        -0.15        0.06    -0.42
> 
> 
> 
> 
> cleanEx(); nameEx("bifactor")
> ### * bifactor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bifactor
> ### Title: Seven data sets showing a bifactor solution.
> ### Aliases: bifactor Bechtoldt.1 Bechtoldt.2 Bechtoldt Holzinger
> ###   Holzinger.9 Reise Thurstone Thurstone.33
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(bifactor)
> holz <- omega(Holzinger,4, title = "14 ability tests from Holzinger-Swineford")
Loading required package: GPArotation
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
> bf <- omega(Reise,5,title="16 health items from Reise") 
> omega(Reise,5,labels=colnames(Reise),title="16 health items from Reise")
16 health items from Reise 
Call: omega(m = Reise, nfactors = 5, title = "16 health items from Reise", 
    labels = colnames(Reise))
Alpha:                 0.91 
G.6:                   0.92 
Omega Hierarchical:    0.82 
Omega H asymptotic:    0.88 
Omega Total            0.93 

Schmid Leiman Factor loadings greater than  0.2 
             g  F1*  F2*  F3*  F4*  F5*   h2   u2
phone     0.63                0.27      0.48 0.52
routine   0.56                0.32      0.42 0.58
illness   0.55                0.25      0.38 0.62
listen    0.77 0.38                     0.74 0.26
explain   0.63 0.39                     0.55 0.45
respect   0.76 0.47                     0.80 0.20
time      0.73 0.37                     0.67 0.33
courtesy  0.68                     0.31 0.59 0.41
helpful   0.77                     0.63 1.00     
happy     0.52      0.22                0.36 0.64
referral  0.44           0.31           0.32 0.68
necessary 0.62           0.54           0.66 0.34
delay     0.52           0.44           0.54 0.46
problem   0.40      0.35                0.29 0.71
help      0.42      0.44                0.38 0.62
paperwork 0.40      0.38                0.30 0.70

With eigenvalues of:
   g  F1*  F2*  F3*  F4*  F5* 
5.79 0.68 0.56 0.61 0.26 0.50 

general/max  8.54   max/min =   2.61
The degrees of freedom for the model is 50  and the fit was  0.02 

Measures of factor score adequacy              g F1* F2* F3* F4* F5*
 Correlation of scores with factors            0.92 0.71 0.61 0.72 0.48 0.89
Multiple R square of scores with factors       0.85 0.5 0.37 0.51 0.23 0.79
Minimum correlation of factor score estimates  0.7 0 -0.26 0.02 -0.54 0.59 
> thur.bf <- omega(Thurstone,title="9 variables from Thurstone")
> 
> 
> 
> 
> cleanEx(); nameEx("circ.tests")
> ### * circ.tests
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: circ.tests
> ### Title: Apply four tests of circumplex versus simple structure
> ### Aliases: circ.tests
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> circ.data <- circ.sim(24,500)
> circ.fa <- factor.pa(circ.data,2)
> #plot(circ.fa$loadings)
> ct <- circ.tests(circ.fa)
> #compare with non-circumplex data
> simp.data <- item.sim(24,500)
> simp.fa <- factor.pa(simp.data,2)
> #plot(simp.fa$loadings)
> st <- circ.tests(simp.fa)
> print(rbind(ct,st),digits=2)
   gaps   fisher RT   VT  
ct 0.0068 0.09   0.35 0.35
st 0.28   0.11   0.89 0.87
> 
> 
> 
> 
> cleanEx(); nameEx("cities")
> ### * cities
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cities
> ### Title: Distances between 11 US cities
> ### Aliases: cities city.location
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> data(cities)
> city.location[,1] <- -city.location[,1]
> if(require(maps)) {map("usa")
+ title("MultiDimensional Scaling of US cities")
+ points(city.location)} else {plot(city.location, xlab="Dimension 1", ylab="Dimension 2",main ="multidimensional scaling of US cities")}
Loading required package: maps
> city.loc <- cmdscale(cities, k=2) #ask for a 2 dimensional solution  round(city.loc,0) 
> city.loc <- -city.loc 
> city.loc <- rescale(city.loc,mean(city.location),sd(city.location))
> points(city.loc,type="n") 
> text(city.loc,labels=names(cities))
> 
> 
> 
> 
> cleanEx(); nameEx("cluster.cor")
> ### * cluster.cor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cluster.cor
> ### Title: Find correlations of composite variables from a larger matrix
> ### Aliases: cluster.cor
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D data(attitude)
> ##D keys <- matrix(c(1,1,1,0,0,0,0,
> ##D                  0,0,0,1,1,1,1),ncol=2)
> ##D colnames(keys) <- c("first","second")
> ##D r.mat <- cor(attitude)
> ##D cluster.cor(keys,r.mat)
> ## End(Not run)
> #$cor
> #       first second
> #first    1.0    0.6
> #second   0.6    1.0
> #
> #$sd
> # first second 
> #  2.57   3.01 
> #
> #$corrected
> #       first second
> #first   0.82   0.77
> #second  0.60   0.74
> #
> #$size
> # first second 
> #     3      4 
> 
> 
> 
> 
> cleanEx(); nameEx("cluster.fit")
> ### * cluster.fit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cluster.fit
> ### Title: cluster Fit: fit of the cluster model to a correlation matrix
> ### Aliases: cluster.fit
> ### Keywords: multivariate cluster
> 
> ### ** Examples
> 
>  r.mat<- Harman74.cor$cov
>  iq.clus <- ICLUST(r.mat,nclusters =2)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
>  fit <- cluster.fit(r.mat,iq.clus$loadings,iq.clus$clusters)
>  fit
$clusterfit
[1] 0.4000123

$structurefit
[1] 0.6466006

$patternfit
[1] 0.8390722

>  
> 
> 
> 
> 
> cleanEx(); nameEx("cluster.loadings")
> ### * cluster.loadings
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cluster.loadings
> ### Title: Find item by cluster correlations, corrected for overlap and
> ###   reliability
> ### Aliases: cluster.loadings
> ### Keywords: multivariate cluster
> 
> ### ** Examples
> 
> 
>  r.mat<- Harman74.cor$cov
>  clusters <- matrix(c(1,1,1,rep(0,24),1,1,1,1,rep(0,17)),ncol=2)
>  cluster.loadings(clusters,r.mat)
Call: cluster.loadings(keys = clusters, r.mat = r.mat)

(Standardized) Alpha:
[1] 0.61 0.79

(Standardized) G6*:
[1] 0.66 0.84

Average item correlation:
[1] 0.35 0.48

Number of items:
[1] 3 4

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
     [,1] [,2]
[1,] 0.61 0.69
[2,] 0.48 0.79

Item by scale intercorrelations
 corrected for item overlap and scale reliability
                        [,1] [,2]
VisualPerception       0.675 0.50
Cubes                  0.512 0.32
PaperFormBoard         0.636 0.36
Flags                  0.550 0.45
GeneralInformation     0.468 0.76
PargraphComprehension  0.459 0.82
SentenceCompletion     0.375 0.84
WordClassification     0.478 0.74
WordMeaning            0.387 0.85
Addition               0.054 0.30
Code                   0.301 0.36
CountingDots           0.328 0.23
StraightCurvedCapitals 0.575 0.46
WordRecognition        0.222 0.30
NumberRecognition      0.238 0.27
FigureRecognition      0.520 0.34
ObjectNumber           0.196 0.31
NumberFigure           0.457 0.29
FigureWord             0.380 0.28
Deduction              0.523 0.57
NumericalPuzzles       0.461 0.43
ProblemReasoning       0.491 0.56
SeriesCompletion       0.661 0.56
ArithmeticProblems     0.382 0.54
> 
>  
> 
> 
> 
> 
> cleanEx(); nameEx("cluster.plot")
> ### * cluster.plot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cluster.plot
> ### Title: Plot factor/cluster loadings and assign items to clusters by
> ###   their highest loading.
> ### Aliases: cluster.plot factor.plot
> ### Keywords: multivariate hplot cluster
> 
> ### ** Examples
> 
> circ.data <- circ.sim(24,500)
> circ.fa <- factor.pa(circ.data,2)
> cluster.plot(circ.fa,cut=.5)
> 
> 
> 
> cleanEx(); nameEx("cluster2keys")
> ### * cluster2keys
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cluster2keys
> ### Title: Convert a cluster vector (from e.g., kmeans) to a keys matrix
> ###   suitable for scoring item clusters.
> ### Aliases: cluster2keys
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> test.data <- Harman74.cor$cov
> kc <- kmeans(test.data,4)
> keys <- cluster2keys(kc)
> keys  #these match those found by ICLUST
      [,1] [,2] [,3] [,4]
 [1,]    1    0    0    0
 [2,]    1    0    0    0
 [3,]    1    0    0    0
 [4,]    1    0    0    0
 [5,]    0    1    0    0
 [6,]    0    1    0    0
 [7,]    0    1    0    0
 [8,]    0    1    0    0
 [9,]    0    1    0    0
[10,]    0    0    1    0
[11,]    0    0    1    0
[12,]    0    0    1    0
[13,]    0    0    1    0
[14,]    0    0    0    1
[15,]    0    0    0    1
[16,]    0    0    0    1
[17,]    0    0    0    1
[18,]    0    0    0    1
[19,]    0    0    0    1
[20,]    1    0    0    0
[21,]    0    0    1    0
[22,]    1    0    0    0
[23,]    1    0    0    0
[24,]    0    0    1    0
> cluster.cor(keys,test.data)
Call: cluster.cor(keys = keys, r.mat = test.data)

(Standardized) Alpha:
[1] 0.80 0.90 0.83 0.74

Average item correlation:
[1] 0.36 0.64 0.45 0.32

Number of items:
[1] 7 5 6 6

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
     [,1] [,2] [,3] [,4]
[1,] 0.80 0.72 0.63 0.69
[2,] 0.61 0.90 0.56 0.52
[3,] 0.51 0.48 0.83 0.67
[4,] 0.53 0.43 0.53 0.74
> 
> 
> 
> cleanEx(); nameEx("comorbidity")
> ### * comorbidity
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: comorbidity
> ### Title: Convert base rates of two diagnoses and their comorbidity into
> ###   phi, Yule, and tetrachorics
> ### Aliases: comorbidity
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> if(require(polycor)) {comorbidity(.2,.15,.1,c("Anxiety","Depression")) }
Loading required package: polycor
Loading required package: mvtnorm
Loading required package: sfsmisc
$twobytwo
            Anxiety -Anxiety
Depression      0.1     0.05
-Depression     0.1     0.75

$phi
[1] 0.49

$Yule
[1] 0.875

$tetra
[1] 0.7518376

> 
> 
> 
> cleanEx(); nameEx("cor.plot")
> ### * cor.plot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cor.plot
> ### Title: Create an image plot for a correlation or factor matrix
> ### Aliases: cor.plot
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> data(bifactor)
> cor.plot(Thurstone,TRUE, main="9 cognitive variables from Thurstone")
> simp <- sim.circ(24)
> cor.plot(cor(simp),colors=TRUE,zlim=c(-1,1),main="24 variables in a circumplex")
> 
> 
> 
> 
> cleanEx(); nameEx("corr.test")
> ### * corr.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: corr.test
> ### Title: Find the correlations, sample sizes, and probability values
> ###   between elements of a matrix or data.frame.
> ### Aliases: corr.test
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> data(sat.act)
> corr.test(sat.act)
Call:corr.test(x = sat.act)
Correlation matrix 
          gender education   age   ACT  SATV  SATQ
gender      1.00      0.09 -0.02 -0.04 -0.02 -0.17
education   0.09      1.00  0.55  0.15  0.05  0.03
age        -0.02      0.55  1.00  0.11 -0.04 -0.03
ACT        -0.04      0.15  0.11  1.00  0.56  0.59
SATV       -0.02      0.05 -0.04  0.56  1.00  0.64
SATQ       -0.17      0.03 -0.03  0.59  0.64  1.00
Sample Size 
          gender education age ACT SATV SATQ
gender       700       700 700 700  700  687
education    700       700 700 700  700  687
age          700       700 700 700  700  687
ACT          700       700 700 700  700  687
SATV         700       700 700 700  700  687
SATQ         687       687 687 687  687  687
Probability value 
          gender education  age  ACT SATV SATQ
gender      0.00      0.02 0.58 0.33 0.62 0.00
education   0.02      0.00 0.00 0.00 0.22 0.36
age         0.58      0.00 0.00 0.00 0.26 0.37
ACT         0.33      0.00 0.00 0.00 0.00 0.00
SATV        0.62      0.22 0.26 0.00 0.00 0.00
SATQ        0.00      0.36 0.37 0.00 0.00 0.00
> 
> 
> 
> cleanEx(); nameEx("correct.cor")
> ### * correct.cor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: correct.cor
> ### Title: Find dis-attenuated correlations given correlations and
> ###   reliabilities
> ### Aliases: correct.cor
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> 
> # attitude from the datasets package
> #example 1 is a rather clunky way of doing things
> 
> a1 <- attitude[,c(1:3)]
> a2 <- attitude[,c(4:7)]
> x1 <- rowSums(a1)  #find the sum of the first 3 attitudes
> x2 <- rowSums(a2)   #find the sum of the last 4 attitudes
> alpha1 <- alpha(a1)
> alpha2 <- alpha(a2)
> x <- matrix(c(x1,x2),ncol=2)
> x.cor <- cor(x)
> alpha <- c(alpha1$total$raw_alpha,alpha2$total$raw_alpha)
> round(correct.cor(x.cor,alpha),2)
     [,1] [,2]
[1,] 0.82 0.78
[2,] 0.61 0.75
> #
> #much better - although uses standardized alpha 
> clusters <- matrix(c(rep(1,3),rep(0,7),rep(1,4)),ncol=2)
> cluster.loadings(clusters,cor(attitude))
Call: cluster.loadings(keys = clusters, r.mat = cor(attitude))

(Standardized) Alpha:
[1] 0.82 0.74

(Standardized) G6*:
[1] 0.83 0.78

Average item correlation:
[1] 0.60 0.42

Number of items:
[1] 3 4

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
     [,1] [,2]
[1,] 0.82 0.77
[2,] 0.60 0.74

Item by scale intercorrelations
 corrected for item overlap and scale reliability
           [,1] [,2]
rating     0.85 0.57
complaints 0.92 0.63
privileges 0.58 0.54
learning   0.73 0.72
raises     0.73 0.85
critical   0.21 0.36
advance    0.31 0.72
> # or 
> clusters <- matrix(c(rep(1,3),rep(0,7),rep(1,4)),ncol=2)
> cluster.cor(clusters,cor(attitude))
Call: cluster.cor(keys = clusters, r.mat = cor(attitude))

(Standardized) Alpha:
[1] 0.82 0.74

Average item correlation:
[1] 0.60 0.42

Number of items:
[1] 3 4

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
     [,1] [,2]
[1,] 0.82 0.77
[2,] 0.60 0.74
> #
> #best
> scores <- score.items(matrix(c(rep(1,3),rep(0,7),rep(1,4)),ncol=2),attitude)
> scores$corrected
     [,1] [,2]
[1,] 0.82 0.78
[2,] 0.61 0.75
> 
> 
> 
> 
> cleanEx(); nameEx("cortest.bartlett")
> ### * cortest.bartlett
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cortest.bartlett
> ### Title: Bartlett's test that a correlation matrix is an identity matrix
> ### Aliases: cortest.bartlett
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> set.seed(42)   
> x <- matrix(rnorm(1000),ncol=10)
> r <- cor(x)
> cortest.bartlett(r)      #random data don't differ from an identity matrix
Warning in cortest.bartlett(r) : n not specified, 100 used
$chisq
[1] 33.78894

$p.value
[1] 0.8897656

$df
[1] 45

> data(bfi)
> cortest.bartlett(bfi)    #not an identity matrix
R was not square, finding R from data
$chisq
[1] 7459.058

$p.value
[1] 0

$df
[1] 300

> 
> 
> 
> cleanEx(); nameEx("cortest.mat")
> ### * cortest.mat
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cortest.mat
> ### Title: Chi square tests of whether a single matrix is an identity
> ###   matrix, or a pair of matrices are equal.
> ### Aliases: cortest.normal cortest.mat cortest.jennrich cortest
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> x <- matrix(rnorm(1000),ncol=10)
> y <- matrix(rnorm(500),ncol=10)
> cortest.normal(x)  #just test if this matrix is an identity
R1 was not square, finding R from data
Tests of correlation matrices 
Call:cortest.normal(R1 = x)
 Chi Square value 40.03  with df =  45   with probability < 0.68> cortest.normal(x,y) #do these two matrices differ?
R1 was not square, finding R from data
R2 was not square, finding R from data
Tests of correlation matrices 
Call:cortest.normal(R1 = x, R2 = y)
 Chi Square value 48.29  with df =  45   with probability < 0.34> cortest.mat(x)
Warning in cortest.mat(x) :
  R1 matrix was not square, correlations found
Bartlett's test of is R = I
Tests of correlation matrices 
Call:cortest.mat(R1 = x)
 Chi Square value 41.05  with df =  45   with probability < 0.64> cortest.mat(x,y)   #twice the degrees of freedom as the Jennrich
Warning in cortest.mat(x, y) :
  R1 matrix was not square, correlations found
Warning in cortest.mat(x, y) :
  R2 matrix was not square, correlations found
Tests of correlation matrices 
Call:cortest.mat(R1 = x, R2 = y)
 Chi Square value 136.83  with df =  90   with probability < 0.0011> cortest.jennrich(x,y) #
Warning in cortest.jennrich(x, y) :
  R1 matrix was not square, correlations found
Warning in cortest.jennrich(x, y) :
  R2 matrix was not square, correlations found
$chi2
[1] 57.1203

$prob
[1] 0.1061990

> 
> 
> 
> 
> cleanEx(); nameEx("cosinor")
> ### * cosinor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cosinor
> ### Title: Functions for analysis of circadian or diurnal data
> ### Aliases: cosinor circadian.mean circadian.cor circadian.linear.cor
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> time <- seq(1:24)
> pure <- matrix(time,24,18)
> pure <- cos((pure + col(pure))*pi/12)
> matplot(pure,type="l") 
> p <- cosinor(time,pure)
> set.seed(42)
> noisey <- pure + rnorm(24*18)
> n <- cosinor(time,noisey) 
> small.pure <- pure[c(6:18),]
> small.noisey <- noisey[c(6:18),]
> sp <- cosinor(time[c(6:18)],small.pure)
> spo <- cosinor(time[c(6:18)],small.pure,opti=TRUE)
> sn <- cosinor(time[c(6:18)],small.noisey)
> sno <- cosinor(time[c(6:18)],small.noisey,opti=TRUE)
> sum.df <- data.frame(pure=p,noisey = n, small=sp,small.noise = sn, small.opt=spo,small.noise.opt=sno)
> round(sum.df,2)
   pure.phase pure.fit noisey.phase noisey.fit small.phase small.fit
1          23        1         2.48       0.50       22.03         1
2          22        1        23.30       0.66       20.62         1
3          21        1        21.03       0.70       19.69         1
4          20        1        19.60       0.72       19.02         1
5          19        1        18.83       0.43       18.48         1
6          18        1        17.13       0.83       18.00         1
7          17        1        18.14       0.69       17.52         1
8          16        1        15.15       0.59       16.98         1
9          15        1        14.13       0.43       16.31         1
10         14        1        14.43       0.70       15.38         1
11         13        1        11.34       0.47       13.97         1
12         12        1        12.72       0.57       12.00         1
13         11        1        11.34       0.72       10.03         1
14         10        1        10.82       0.74        8.62         1
15          9        1        10.52       0.72        7.69         1
16          8        1         8.63       0.48        7.02         1
17          7        1         7.50       0.81        6.48         1
18          6        1         4.78       0.31        6.00         1
   small.noise.phase small.noise.fit small.opt.phase small.opt.fit
1               7.33            0.41              23             1
2              23.95            0.54              22             1
3              21.04            0.71              21             1
4              19.52            0.76              20             1
5              18.92            0.60              19             1
6              17.70            0.84              18             1
7              17.16            0.72              17             1
8              16.68            0.40              16             1
9              12.66            0.16              15             1
10             14.59            0.75              14             1
11              3.80            0.27              13             1
12             11.84            0.54              12             1
13              6.32            0.37              11             1
14              9.96            0.63              10             1
15              8.45            0.57               9             1
16              7.91            0.68               8             1
17              6.49            0.72               7             1
18              7.33            0.39               6             1
   small.noise.opt.phase small.noise.opt.fit
1                   8.50                0.41
2                  23.98                0.54
3                  22.34                0.71
4                  20.78                0.76
5                  19.84                0.60
6                  17.36                0.84
7                  16.31                0.72
8                  15.51                0.40
9                  12.32                0.16
10                 13.39                0.75
11                  2.40                0.27
12                 11.92                0.54
13                  6.67                0.37
14                 10.96                0.63
15                  9.85                0.57
16                  9.27                0.68
17                  7.02                0.72
18                  8.51                0.39
> round(circadian.cor(sum.df[,c(1,3,5,7,9,11)]),2)  #compare alternatives 
                      pure.phase noisey.phase small.phase small.noise.phase
pure.phase                  1.00         0.93        0.98              0.68
noisey.phase                0.93         1.00        0.87              0.87
small.phase                 0.98         0.87        1.00              0.62
small.noise.phase           0.68         0.87        0.62              1.00
small.opt.phase             1.00         0.93        0.98              0.68
small.noise.opt.phase       0.68         0.86        0.62              0.99
                      small.opt.phase small.noise.opt.phase
pure.phase                       1.00                  0.68
noisey.phase                     0.93                  0.86
small.phase                      0.98                  0.62
small.noise.phase                0.68                  0.99
small.opt.phase                  1.00                  0.68
small.noise.opt.phase            0.68                  1.00
> round(cor(sum.df[,c(2,4,6,8,10,12)]),2) 
                    pure.fit noisey.fit small.fit small.noise.fit small.opt.fit
pure.fit                1.00      -0.22      0.00            0.15          0.09
noisey.fit             -0.22       1.00      0.18            0.65          0.06
small.fit               0.00       0.18      1.00            0.13         -0.14
small.noise.fit         0.15       0.65      0.13            1.00          0.48
small.opt.fit           0.09       0.06     -0.14            0.48          1.00
small.noise.opt.fit     0.15       0.65      0.13            1.00          0.48
                    small.noise.opt.fit
pure.fit                           0.15
noisey.fit                         0.65
small.fit                          0.13
small.noise.fit                    1.00
small.opt.fit                      0.48
small.noise.opt.fit                1.00
> 
> 
> 
> cleanEx(); nameEx("count.pairwise")
> ### * count.pairwise
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: count.pairwise
> ### Title: Count number of pairwise cases for a data set with missing (NA)
> ###   data.
> ### Aliases: count.pairwise
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> ## Not run: 
> ##D x <- matrix(rnorm(1000),ncol=6)
> ##D y <- matrix(rnorm(500),ncol=3)
> ##D x[x < 0] <- NA
> ##D y[y> 1] <- NA
> ##D 
> ##D count.pairwise(x)
> ##D count.pairwise(y)
> ##D count.pairwise(x,y)
> ## End(Not run)
> 
>     
> 
> 
> 
> cleanEx(); nameEx("cta")
> ### * cta
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cta
> ### Title: Simulate the C(ues) T(endency) A(ction) model of motivation
> ### Aliases: cta
> ### Keywords: models
> 
> ### ** Examples
> 
> #not run 
> #cta()   #default values, running over time 
> #cta(type="state") #default values, in a state space  of tendency 1 versus tendency 2
> 
> 
> 
> cleanEx(); nameEx("cubits")
> ### * cubits
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cubits
> ### Title: Galton's example of the relationship between height and 'cubit'
> ###   or forearm length
> ### Aliases: cubits
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(cubits)
> cubits
   16.5 16.75 17.25 17.75 18.25 18.75 19.25 19.75
71    0     0     0     1     3     4    15     7
70    0     0     0     1     5    13    11     0
69    0     1     1     2    25    15     6     0
68    0     1     3     7    14     7     4     2
67    0     1     7    15    28     8     2     0
66    0     1     7    18    15     6     0     0
65    0     4    10    12     8     2     0     0
64    0     5    11     2     3     0     0     0
63    9    12    10     3     1     0     0     0
> heights <- table2df(cubits,labs <- c("height","cubit"))
> ellipses(heights,n=1,main="Galton's co-relation data set")
> ellipses(jitter(heights$cubit,3),jitter(heights$height,3),pch=".",main="Galton's co-relation data set") #add in some noise to see the points
> 
> 
> 
> 
> cleanEx(); nameEx("describe")
> ### * describe
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: describe
> ### Title: Basic descriptive statistics useful for psychometrics
> ### Aliases: describe
> ### Keywords: multivariate models univar
> 
> ### ** Examples
> 
> data(sat.act)
> describe(sat.act)
          var   n   mean     sd median trimmed    mad min max range  skew
gender      1 700   1.65   0.48      2    1.68   0.00   1   2     1 -0.61
education   2 700   3.16   1.43      3    3.31   1.48   0   5     5 -0.68
age         3 700  25.59   9.50     22   23.86   5.93  13  65    52  1.64
ACT         4 700  28.55   4.82     29   28.84   4.45   3  36    33 -0.66
SATV        5 700 612.23 112.90    620  619.45 118.61 200 800   600 -0.64
SATQ        6 687 610.22 115.64    620  617.25 118.61 200 800   600 -0.59
          kurtosis   se
gender       -1.62 0.02
education    -0.07 0.05
age           2.42 0.36
ACT           0.53 0.18
SATV          0.33 4.27
SATQ         -0.02 4.41
> 
> describe(sat.act,skew=FALSE)   
          var   n   mean     sd median trimmed    mad min max range   se
gender      1 700   1.65   0.48      2    1.68   0.00   1   2     1 0.02
education   2 700   3.16   1.43      3    3.31   1.48   0   5     5 0.05
age         3 700  25.59   9.50     22   23.86   5.93  13  65    52 0.36
ACT         4 700  28.55   4.82     29   28.84   4.45   3  36    33 0.18
SATV        5 700 612.23 112.90    620  619.45 118.61 200 800   600 4.27
SATQ        6 687 610.22 115.64    620  617.25 118.61 200 800   600 4.41
> 
> 
> 
> 
> 
> 
> 
> cleanEx(); nameEx("describe.by")
> ### * describe.by
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: describe.by
> ### Title: Basic summary statistics by group
> ### Aliases: describe.by
> ### Keywords: models univar
> 
> ### ** Examples
> 
> data(sat.act)
> describe.by(sat.act,sat.act$gender) #just one grouping variable 
$`1`
          var   n   mean     sd median trimmed    mad min max range  skew
gender      1 247   1.00   0.00      1    1.00   0.00   1   1     0   NaN
education   2 247   3.00   1.54      3    3.12   1.48   0   5     5 -0.54
age         3 247  25.86   9.74     22   24.23   5.93  14  58    44  1.43
ACT         4 247  28.79   5.06     30   29.23   4.45   3  36    33 -1.06
SATV        5 247 615.11 114.16    630  622.07 118.61 200 800   600 -0.63
SATQ        6 245 635.87 116.02    660  645.53  94.89 300 800   500 -0.72
          kurtosis   se
gender         NaN 0.00
education    -0.60 0.10
age           1.43 0.62
ACT           1.89 0.32
SATV          0.13 7.26
SATQ         -0.12 7.41

$`2`
          var   n   mean     sd median trimmed    mad min max range  skew
gender      1 453   2.00   0.00      2    2.00   0.00   2   2     0   NaN
education   2 453   3.26   1.35      3    3.40   1.48   0   5     5 -0.74
age         3 453  25.45   9.37     22   23.70   5.93  13  65    52  1.77
ACT         4 453  28.42   4.69     29   28.63   4.45  15  36    21 -0.39
SATV        5 453 610.66 112.31    620  617.91 103.78 200 800   600 -0.65
SATQ        6 442 596.00 113.07    600  602.21 133.43 200 800   600 -0.58
          kurtosis   se
gender         NaN 0.00
education     0.27 0.06
age           3.03 0.44
ACT          -0.42 0.22
SATV          0.42 5.28
SATQ          0.13 5.38

> #describe.by(sat.act,list(sat.act$gender,sat.act$education))  #two grouping variables
> #des.mat <- describe.by(sat.act$age,sat.act$education,mat=TRUE) #matrix output 
> #des.mat <- describe.by(sat.act$age,list(sat.act$education,sat.act$gender),mat=TRUE) 
> 
> 
> 
> 
> cleanEx(); nameEx("eigen.loadings")
> ### * eigen.loadings
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: eigen.loadings
> ### Title: Convert eigen vectors and eigen values to the more normal (for
> ###   psychologists) component loadings
> ### Aliases: eigen.loadings
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> x <- eigen(Harman74.cor$cov)
> x$vectors[1:8,1:4]  #as they appear from eigen
           [,1]         [,2]       [,3]        [,4]
[1,] -0.2158754 -0.003763753 -0.3287459  0.16684940
[2,] -0.1401069 -0.054849237 -0.3075102  0.16446394
[3,] -0.1558742 -0.131808462 -0.3661446  0.08629700
[4,] -0.1789834 -0.122936419 -0.2572906  0.17620811
[5,] -0.2436443 -0.221950283  0.2577381  0.04309103
[6,] -0.2420622 -0.288461492  0.2037767 -0.06579018
[7,] -0.2372958 -0.293489427  0.2731887  0.05917239
[8,] -0.2433556 -0.167723001  0.1103619  0.09493344
> y <- princomp(covmat=Harman74.cor$cov) 
> y$loadings[1:8,1:4] #as they appear from princomp
                          Comp.1       Comp.2     Comp.3      Comp.4
VisualPerception      -0.2158754 -0.003763753 -0.3287459  0.16684940
Cubes                 -0.1401069 -0.054849237 -0.3075102  0.16446394
PaperFormBoard        -0.1558742 -0.131808462 -0.3661446  0.08629700
Flags                 -0.1789834 -0.122936419 -0.2572906  0.17620811
GeneralInformation    -0.2436443 -0.221950283  0.2577381  0.04309103
PargraphComprehension -0.2420622 -0.288461492  0.2037767 -0.06579018
SentenceCompletion    -0.2372958 -0.293489427  0.2731887  0.05917239
WordClassification    -0.2433556 -0.167723001  0.1103619  0.09493344
> eigen.loadings(x)[1:8,1:4] # rescaled by the eigen values
           [,1]         [,2]       [,3]        [,4]
[1,] -0.6157349 -0.005449052 -0.4276989  0.20447285
[2,] -0.3996226 -0.079409132 -0.4000712  0.20154949
[3,] -0.4445954 -0.190828463 -0.4763546  0.10575641
[4,] -0.5105089 -0.177983777 -0.3347354  0.21594189
[5,] -0.6949395 -0.321333174  0.3353177  0.05280778
[6,] -0.6904266 -0.417626171  0.2651138 -0.08062544
[7,] -0.6768317 -0.424905469  0.3554189  0.07251538
[8,] -0.6941158 -0.242824490  0.1435811  0.11634031
> 
> 
> 
> cleanEx(); nameEx("ellipses")
> ### * ellipses
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ellipses
> ### Title: Plot data and 1 and 2 sigma correlation ellipses
> ### Aliases: ellipses
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> data(galton)
> ellipses(galton,lm=TRUE)
> ellipses(galton$parent,galton$child,xlab="Mid Parent Height",ylab="Child Height") #input are two vectors
> data(sat.act)
> ellipses(sat.act)  #shows the pairs.panels ellipses
> 
> 
> 
> cleanEx(); nameEx("epi.bfi")
> ### * epi.bfi
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: epi.bfi
> ### Title: 13 personality scales from the Eysenck Personality Inventory and
> ###   Big 5 inventory
> ### Aliases: epi.bfi
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(epi.bfi)
> pairs.panels(epi.bfi[,1:5])
> describe(epi.bfi)
         var   n   mean    sd median trimmed   mad min max range  skew kurtosis
epiE       1 231  13.33  4.14     14   13.49  4.45   1  22    21 -0.33    -0.06
epiS       2 231   7.58  2.69      8    7.77  2.97   0  13    13 -0.57    -0.02
epiImp     3 231   4.37  1.88      4    4.36  1.48   0   9     9  0.06    -0.62
epilie     4 231   2.38  1.50      2    2.27  1.48   0   7     7  0.66     0.24
epiNeur    5 231  10.41  4.90     10   10.39  4.45   0  23    23  0.06    -0.50
bfagree    6 231 125.00 18.14    126  125.26 17.79  74 167    93 -0.21    -0.27
bfcon      7 231 113.25 21.88    114  113.42 22.24  53 178   125 -0.02     0.23
bfext      8 231 102.18 26.45    104  102.99 22.24   8 168   160 -0.41     0.51
bfneur     9 231  87.97 23.34     90   87.70 23.72  34 152   118  0.07    -0.55
bfopen    10 231 123.43 20.51    125  123.78 20.76  73 173   100 -0.16    -0.16
bdi       11 231   6.78  5.78      6    5.97  4.45   0  27    27  1.29     1.50
traitanx  12 231  39.01  9.52     38   38.36  8.90  22  71    49  0.67     0.47
stateanx  13 231  39.85 11.48     38   38.92 10.38  21  79    58  0.72    -0.01
           se
epiE     0.27
epiS     0.18
epiImp   0.12
epilie   0.10
epiNeur  0.32
bfagree  1.19
bfcon    1.44
bfext    1.74
bfneur   1.54
bfopen   1.35
bdi      0.38
traitanx 0.63
stateanx 0.76
> 
> 
> 
> cleanEx(); nameEx("error.bars")
> ### * error.bars
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: error.bars
> ### Title: Plot means and confidence intervals
> ### Aliases: error.bars
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> x <- replicate(20,rnorm(50))
> boxplot(x,notch=TRUE,main="Notched boxplot with error bars")
> error.bars(x,add=TRUE)
> abline(h=0)
> 
> error.bars(attitude,alpha=.5,main="50 percent confidence limits") #another example
> error.bars(attitude,bar=TRUE)  #show the use of bar graphs
> 
> #combine with a strip chart
> stripchart(attitude,vertical=TRUE,method="jitter",main="Stripchart with 95 percent confidence limits")
> error.bars(attitude,add=TRUE,arrow.len=.2)
> 
> 
> 
> 
> cleanEx(); nameEx("error.bars.by")
> ### * error.bars.by
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: error.bars.by
> ### Title: Plot means and confidence intervals for multiple groups
> ### Aliases: error.bars.by
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> 
> data(sat.act)
> error.bars.by(sat.act[1:4],sat.act$gender)
> error.bars.by(sat.act[5:6],sat.act$gender,bars=TRUE,labels=c("male","female"),main="SAT V and SAT Q by gender")  #draw a barplot
> 
> error.bars.by(sat.act[5:6],sat.act$education,bars=TRUE,xlab="Education",main="95 percent confidence limits of Sat V and Sat Q") 
> 
> error.bars.by(sat.act[5:6],sat.act$education,TRUE, xlab="Education")  #plot SAT V and SAT Q by education
> 
> 
> 
> 
> cleanEx(); nameEx("error.crosses")
> ### * error.crosses
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: error.crosses
> ### Title: Plot x and y error bars
> ### Aliases: error.crosses
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> 
> desc <- describe(attitude)
> x <- desc[1,]
> y <- desc[2,]
> plot(x$mean,y$mean,xlab=rownames(x),ylab=rownames(y))   #in graphics window
> error.crosses(x,y)    #in graphics window
> #now for a bit more complicated plotting 
> desc <- describe.by(attitude,(attitude[,7]>41)) #select a high and low group
> g1 <- desc$'FALSE'
> g2 <- desc$'TRUE'
> plot(g1$mean,g2$mean,xlab = "Low Advance",ylab="High Advance",xlim=c(30,80),ylim=c(50,80))
> error.crosses(g1,g2,labels=rownames(g1),pos=rep(1,7))
> title("Attitudes grouped by high and low scores on Advance")
> 
> 
> 
> 
> cleanEx(); nameEx("fa")
> ### * fa
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fa
> ### Title: Factor analysis by Principal Axis, MinRes (minimum residual),
> ###   Weighted Least Squares or Maximum Likelihood
> ### Aliases: fa factor.pa factor.minres factor.wls
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> #using the Harman 24 mental tests, compare a principal factor with a principal components solution
> pc <- principal(Harman74.cor$cov,4,rotate="varimax")
> pa <- fa(Harman74.cor$cov,4,fm="pa" ,rotate="varimax")
> uls <- fa(Harman74.cor$cov,4,rotate="varimax")
> wls <- fa(Harman74.cor$cov,4,fm="wls")
> 
> #to show the loadings sorted by absolute value
> print(uls,sort=TRUE)
Factor Analysis using method =  minres
Call: fa(r = Harman74.cor$cov, nfactors = 4, rotate = "varimax")
                        V  MR1   MR3   MR2  MR4   h2   u2
SentenceCompletion      7 0.81                  0.72 0.28
WordMeaning             9 0.81                  0.74 0.26
PargraphComprehension   6 0.76                  0.68 0.32
GeneralInformation      5 0.74                  0.64 0.36
WordClassification      8 0.57  0.34            0.52 0.48
VisualPerception        1       0.69            0.56 0.44
PaperFormBoard          3       0.57            0.35 0.65
Flags                   4       0.53            0.35 0.65
SeriesCompletion       23 0.37  0.51            0.51 0.49
Cubes                   2       0.44            0.22 0.78
Deduction              20 0.38  0.41            0.41 0.59
ProblemReasoning       22 0.37  0.40            0.40 0.60
Addition               10             0.82      0.75 0.25
CountingDots           12             0.71      0.56 0.44
StraightCurvedCapitals 13       0.43  0.53      0.51 0.49
Code                   11             0.52 0.38 0.46 0.54
ArithmeticProblems     24 0.37        0.50      0.50 0.50
NumericalPuzzles       21       0.39  0.44      0.42 0.58
ObjectNumber           17                  0.58 0.41 0.59
WordRecognition        14                  0.56 0.36 0.64
FigureRecognition      16       0.41       0.52 0.45 0.55
NumberRecognition      15                  0.52 0.31 0.69
NumberFigure           18             0.34 0.45 0.41 0.59
FigureWord             19                  0.36 0.24 0.76

                MR1  MR3  MR2  MR4
SS loadings    3.64 2.89 2.66 2.28
Proportion Var 0.15 0.12 0.11 0.09
Cumulative Var 0.15 0.27 0.38 0.48

Test of the hypothesis that 4 factors are sufficient.

The degrees of freedom for the model is 186  and the fit was  1.71 

Measures of factor score adequacy              MR1 MR3 MR2 MR4
 Correlation of scores with factors            0.93 0.87 0.91 0.83
Multiple R square of scores with factors       0.87 0.76 0.82 0.68
Minimum correlation of factor score estimates  0.73 0.51 0.65 0.36
Validity of unit weighted factor scores        0.87 0.75 0.8 0.76 
> 
> #then compare with a maximum likelihood solution using factanal
> mle <- factanal(covmat=Harman74.cor$cov,factors=4)
> factor.congruence(list(mle,pa,pc,uls,wls))
        Factor1 Factor2 Factor3 Factor4  PA1  PA3  PA2  PA4  PC1  PC3  PC2  PC4
Factor1    1.00    0.61    0.46    0.56 1.00 0.61 0.46 0.55 1.00 0.54 0.44 0.47
Factor2    0.61    1.00    0.50    0.61 0.61 1.00 0.50 0.60 0.60 0.99 0.49 0.52
Factor3    0.46    0.50    1.00    0.57 0.46 0.50 1.00 0.56 0.45 0.44 1.00 0.48
Factor4    0.56    0.61    0.57    1.00 0.56 0.62 0.58 1.00 0.55 0.55 0.56 0.99
PA1        1.00    0.61    0.46    0.56 1.00 0.61 0.46 0.55 1.00 0.54 0.44 0.47
PA3        0.61    1.00    0.50    0.62 0.61 1.00 0.50 0.61 0.61 0.99 0.50 0.53
PA2        0.46    0.50    1.00    0.58 0.46 0.50 1.00 0.57 0.46 0.44 1.00 0.49
PA4        0.55    0.60    0.56    1.00 0.55 0.61 0.57 1.00 0.54 0.54 0.55 0.99
PC1        1.00    0.60    0.45    0.55 1.00 0.61 0.46 0.54 1.00 0.53 0.43 0.46
PC3        0.54    0.99    0.44    0.55 0.54 0.99 0.44 0.54 0.53 1.00 0.43 0.47
PC2        0.44    0.49    1.00    0.56 0.44 0.50 1.00 0.55 0.43 0.43 1.00 0.47
PC4        0.47    0.52    0.48    0.99 0.47 0.53 0.49 0.99 0.46 0.47 0.47 1.00
MR1        1.00    0.61    0.46    0.56 1.00 0.61 0.46 0.55 1.00 0.54 0.44 0.47
MR3        0.61    1.00    0.50    0.61 0.61 1.00 0.50 0.60 0.60 0.99 0.49 0.53
MR2        0.46    0.50    1.00    0.58 0.46 0.51 1.00 0.57 0.46 0.44 1.00 0.48
MR4        0.55    0.61    0.57    1.00 0.55 0.61 0.58 1.00 0.55 0.55 0.56 0.99
WLS1       1.00    0.61    0.46    0.56 1.00 0.62 0.46 0.55 1.00 0.54 0.44 0.47
WLS3       0.61    1.00    0.50    0.62 0.61 1.00 0.51 0.61 0.61 0.99 0.50 0.53
WLS2       0.47    0.50    1.00    0.58 0.47 0.51 1.00 0.58 0.46 0.44 1.00 0.49
WLS4       0.54    0.59    0.56    1.00 0.54 0.60 0.56 1.00 0.54 0.53 0.55 0.99
         MR1  MR3  MR2  MR4 WLS1 WLS3 WLS2 WLS4
Factor1 1.00 0.61 0.46 0.55 1.00 0.61 0.47 0.54
Factor2 0.61 1.00 0.50 0.61 0.61 1.00 0.50 0.59
Factor3 0.46 0.50 1.00 0.57 0.46 0.50 1.00 0.56
Factor4 0.56 0.61 0.58 1.00 0.56 0.62 0.58 1.00
PA1     1.00 0.61 0.46 0.55 1.00 0.61 0.47 0.54
PA3     0.61 1.00 0.51 0.61 0.62 1.00 0.51 0.60
PA2     0.46 0.50 1.00 0.58 0.46 0.51 1.00 0.56
PA4     0.55 0.60 0.57 1.00 0.55 0.61 0.58 1.00
PC1     1.00 0.60 0.46 0.55 1.00 0.61 0.46 0.54
PC3     0.54 0.99 0.44 0.55 0.54 0.99 0.44 0.53
PC2     0.44 0.49 1.00 0.56 0.44 0.50 1.00 0.55
PC4     0.47 0.53 0.48 0.99 0.47 0.53 0.49 0.99
MR1     1.00 0.61 0.46 0.55 1.00 0.61 0.47 0.54
MR3     0.61 1.00 0.50 0.61 0.61 1.00 0.50 0.60
MR2     0.46 0.50 1.00 0.57 0.46 0.51 1.00 0.56
MR4     0.55 0.61 0.57 1.00 0.56 0.62 0.58 1.00
WLS1    1.00 0.61 0.46 0.56 1.00 0.62 0.47 0.55
WLS3    0.61 1.00 0.51 0.62 0.62 1.00 0.51 0.60
WLS2    0.47 0.50 1.00 0.58 0.47 0.51 1.00 0.57
WLS4    0.54 0.60 0.56 1.00 0.55 0.60 0.57 1.00
> #note that the order of factors and the sign of some of factors differ 
> 
> #finally, compare the unrotated factor, ml, uls, and  wls solutions
> wls <- factor.wls(Harman74.cor$cov,4,rotate="none")
> pa <- factor.pa(Harman74.cor$cov,4,rotate="none")
> mle <-  factanal(factors=4,covmat=Harman74.cor$cov,rotation="none")
> uls <- factor.minres(Harman74.cor$cov,4,rotate="none")
> 
> factor.congruence(list(mle,pa,uls,wls))
        Factor1 Factor2 Factor3 Factor4  PA1   PA2   PA3   PA4   MR1   MR2  MR3
Factor1    1.00    0.11    0.25    0.06 1.00 -0.04 -0.05 -0.01  1.00  0.13 0.24
Factor2    0.11    1.00    0.06    0.07 0.14  0.98 -0.08 -0.08  0.11  1.00 0.01
Factor3    0.25    0.06    1.00    0.01 0.30  0.10  0.95  0.02  0.25  0.11 1.00
Factor4    0.06    0.07    0.01    1.00 0.07  0.13 -0.04  0.99  0.06  0.08 0.02
PA1        1.00    0.14    0.30    0.07 1.00  0.00  0.00  0.00  1.00  0.17 0.28
PA2       -0.04    0.98    0.10    0.13 0.00  1.00  0.00  0.00 -0.04  0.98 0.04
PA3       -0.05   -0.08    0.95   -0.04 0.00  0.00  1.00  0.00 -0.05 -0.03 0.96
PA4       -0.01   -0.08    0.02    0.99 0.00  0.00  0.00  1.00 -0.01 -0.07 0.04
MR1        1.00    0.11    0.25    0.06 1.00 -0.04 -0.05 -0.01  1.00  0.13 0.24
MR2        0.13    1.00    0.11    0.08 0.17  0.98 -0.03 -0.07  0.13  1.00 0.05
MR3        0.24    0.01    1.00    0.02 0.28  0.04  0.96  0.04  0.24  0.05 1.00
MR4        0.05    0.05    0.00    1.00 0.06  0.12 -0.06  0.99  0.05  0.06 0.00
WLS1       1.00    0.14    0.30    0.07 1.00  0.00  0.00  0.00  1.00  0.17 0.29
WLS2      -0.04    0.98    0.09    0.15 0.00  1.00 -0.01  0.01 -0.04  0.98 0.03
WLS3      -0.05   -0.07    0.95   -0.04 0.00  0.01  1.00  0.01 -0.05 -0.03 0.96
WLS4      -0.01   -0.09    0.02    0.98 0.00 -0.01  0.00  1.00 -0.01 -0.08 0.03
          MR4 WLS1  WLS2  WLS3  WLS4
Factor1  0.05 1.00 -0.04 -0.05 -0.01
Factor2  0.05 0.14  0.98 -0.07 -0.09
Factor3  0.00 0.30  0.09  0.95  0.02
Factor4  1.00 0.07  0.15 -0.04  0.98
PA1      0.06 1.00  0.00  0.00  0.00
PA2      0.12 0.00  1.00  0.01 -0.01
PA3     -0.06 0.00 -0.01  1.00  0.00
PA4      0.99 0.00  0.01  0.01  1.00
MR1      0.05 1.00 -0.04 -0.05 -0.01
MR2      0.06 0.17  0.98 -0.03 -0.08
MR3      0.00 0.29  0.03  0.96  0.03
MR4      1.00 0.06  0.13 -0.05  0.99
WLS1     0.06 1.00  0.00  0.00  0.00
WLS2     0.13 0.00  1.00  0.00  0.00
WLS3    -0.05 0.00  0.00  1.00  0.00
WLS4     0.99 0.00  0.00  0.00  1.00
> #note that the order of factors and the sign of some of factors differ 
> 
> 
> #an example of where the ML and PA and MR models differ is found in Thurstone.33.
> #compare the first two factors with the 3 factor solution 
> data(bifactor)
> Thurstone.33 <- as.matrix(Thurstone.33)
> mle2 <- factanal(covmat=Thurstone.33,factors=2,rotation="none")
> mle3 <- factanal(covmat=Thurstone.33,factors=3 ,rotation="none")
> pa2 <- factor.pa(Thurstone.33,2,rotate="none")
> pa3 <- factor.pa(Thurstone.33,3,rotate="none")
> mr2 <- fa(Thurstone.33,2,rotate="none")
> mr3 <- fa(Thurstone.33,3,rotate="none")
> factor.congruence(list(mle2,mle3,pa2,pa3,mr2,mr3))
        Factor1 Factor2 Factor1 Factor2 Factor3  PA1   PA2   PA1   PA2   PA3
Factor1    1.00    0.16    0.98    0.88    0.16 1.00 -0.03  1.00 -0.02  0.00
Factor2    0.16    1.00    0.03    0.51   -0.79 0.19  0.98  0.18  0.97 -0.09
Factor1    0.98    0.03    1.00    0.76    0.19 0.97 -0.14  0.98 -0.17 -0.14
Factor2    0.88    0.51    0.76    1.00   -0.01 0.89  0.32  0.88  0.39  0.26
Factor3    0.16   -0.79    0.19   -0.01    1.00 0.15 -0.87  0.15 -0.78  0.61
PA1        1.00    0.19    0.97    0.89    0.15 1.00  0.00  1.00  0.01  0.01
PA2       -0.03    0.98   -0.14    0.32   -0.87 0.00  1.00 -0.01  0.98 -0.17
PA1        1.00    0.18    0.98    0.88    0.15 1.00 -0.01  1.00  0.00  0.00
PA2       -0.02    0.97   -0.17    0.39   -0.78 0.01  0.98  0.00  1.00  0.00
PA3        0.00   -0.09   -0.14    0.26    0.61 0.01 -0.17  0.00  0.00  1.00
MR1        1.00    0.18    0.97    0.89    0.15 1.00 -0.01  1.00  0.00  0.01
MR2        0.03    0.99   -0.09    0.38   -0.85 0.06  1.00  0.05  0.98 -0.14
MR1        0.98    0.05    1.00    0.79    0.18 0.98 -0.13  0.98 -0.14 -0.11
MR2        0.84    0.56    0.72    1.00   -0.05 0.86  0.38  0.85  0.45  0.28
MR3        0.17   -0.77    0.19    0.01    1.00 0.16 -0.85  0.16 -0.76  0.63
          MR1   MR2   MR1   MR2   MR3
Factor1  1.00  0.03  0.98  0.84  0.17
Factor2  0.18  0.99  0.05  0.56 -0.77
Factor1  0.97 -0.09  1.00  0.72  0.19
Factor2  0.89  0.38  0.79  1.00  0.01
Factor3  0.15 -0.85  0.18 -0.05  1.00
PA1      1.00  0.06  0.98  0.86  0.16
PA2     -0.01  1.00 -0.13  0.38 -0.85
PA1      1.00  0.05  0.98  0.85  0.16
PA2      0.00  0.98 -0.14  0.45 -0.76
PA3      0.01 -0.14 -0.11  0.28  0.63
MR1      1.00  0.05  0.98  0.85  0.16
MR2      0.05  1.00 -0.07  0.44 -0.83
MR1      0.98 -0.07  1.00  0.74  0.19
MR2      0.85  0.44  0.74  1.00 -0.03
MR3      0.16 -0.83  0.19 -0.03  1.00
> 
> 
> 
> cleanEx(); nameEx("fa.graph")
> ### * fa.graph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fa.graph
> ### Title: Graph factor loading matrices
> ### Aliases: fa.graph
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> 
> test.simple <- factor.pa(item.sim(16),2)
> if(require(Rgraphviz)) {fa.graph(test.simple) }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
A graphNEL graph with directed edges
Number of Nodes = 18 
Number of Edges = 16 
> 
> 
> 
> 
> cleanEx(); nameEx("fa.parallel")
> ### * fa.parallel
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fa.parallel
> ### Title: Scree plots of data or correlation matrix compared to random
> ###   ``parallel" matrices
> ### Aliases: fa.parallel
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> 
> test.data <- Harman74.cor$cov 
> fa.parallel(test.data,n.obs=200)
> 
> fa.parallel(attitude) 
> #
> 
> 
> 
> 
> cleanEx(); nameEx("factor.congruence")
> ### * factor.congruence
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor.congruence
> ### Title: Coefficient of factor congruence
> ### Aliases: factor.congruence
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> #fa <- factanal(x,4,covmat=Harman74.cor$cov)
> #pc <- principal(Harman74.cor$cov,4)
> #pa <- factor.pa(Harman74.cov$cor,4)
> #factor.congruence(fa,pc)
> #
> #    Factor1 Factor2 Factor3 Factor4
> #PC1    1.00    0.60    0.45    0.55
> #PC2    0.44    0.49    1.00    0.56
> #PC3    0.54    0.99    0.44    0.55
> #PC4    0.47    0.52    0.48    0.99
> 
> # pa <- factor.pa(Harman74.cor$cov,4)
> # factor.congruence(fa,pa)
> #         PA1  PA3  PA2  PA4
> #Factor1 1.00 0.61 0.46 0.55
> #Factor2 0.61 1.00 0.50 0.60
> #Factor3 0.46 0.50 1.00 0.57
> #Factor4 0.56 0.62 0.58 1.00
> 
> #compare with 
> #round(cor(fa$loading,pc$loading),2)
> 
> 
> 
> 
> cleanEx(); nameEx("factor.fit")
> ### * factor.fit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor.fit
> ### Title: How well does the factor model fit a correlation matrix. Part of
> ###   the VSS package
> ### Aliases: factor.fit
> ### Keywords: models models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D #compare the fit of 4 to 3 factors for the Harman 24 variables
> ##D fa4 <- factanal(x,4,covmat=Harman74.cor$cov)
> ##D round(factor.fit(Harman74.cor$cov,fa4$loading),2)
> ##D #[1] 0.9
> ##D fa3 <- factanal(x,3,covmat=Harman74.cor$cov)
> ##D round(factor.fit(Harman74.cor$cov,fa3$loading),2)
> ##D #[1] 0.88
> ##D 
> ## End(Not run)
> 
> 
> 
> 
> cleanEx(); nameEx("factor.model")
> ### * factor.model
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor.model
> ### Title: Find R = F F' + U2 is the basic factor model
> ### Aliases: factor.model
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> f2 <- matrix(c(.9,.8,.7,rep(0,6),.6,.7,.8),ncol=2)
> mod <- factor.model(f2)
> round(mod,2)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.81 0.72 0.63 0.00 0.00 0.00
[2,] 0.72 0.64 0.56 0.00 0.00 0.00
[3,] 0.63 0.56 0.49 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.36 0.42 0.48
[5,] 0.00 0.00 0.00 0.42 0.49 0.56
[6,] 0.00 0.00 0.00 0.48 0.56 0.64
> 
> 
> 
> cleanEx(); nameEx("factor.residuals")
> ### * factor.residuals
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor.residuals
> ### Title: R* = R- F F'
> ### Aliases: factor.residuals
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> fa2 <- factor.pa(Harman74.cor$cov,2,rotate=TRUE)
>  fa2resid <- factor.residuals(Harman74.cor$cov,fa2)
>  fa2resid[1:4,1:4] #residuals with two factors extracted
                 VisualPerception     Cubes PaperFormBoard     Flags
VisualPerception        0.6534250 0.1017699      0.1656195 0.1920302
Cubes                   0.1017699 0.8632542      0.1637670 0.0522387
PaperFormBoard          0.1656195 0.1637670      0.8230977 0.1004048
Flags                   0.1920302 0.0522387      0.1004048 0.7633040
>  fa4 <- factor.pa(Harman74.cor$cov,4,rotate=TRUE)
>  fa4resid <- factor.residuals(Harman74.cor$cov,fa4)
>  fa4resid[1:4,1:4] #residuals with 4 factors extracted
                 VisualPerception       Cubes PaperFormBoard       Flags
VisualPerception       0.44988061 -0.03531723    -0.01017231  0.04142422
Cubes                 -0.03531723  0.77011860     0.04334670 -0.04976063
PaperFormBoard        -0.01017231  0.04334670     0.66147936 -0.02954846
Flags                  0.04142422 -0.04976063    -0.02954846  0.65010907
> 
> 
> 
> 
> cleanEx(); nameEx("factor.rotate")
> ### * factor.rotate
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor.rotate
> ### Title: ``Hand" rotate a factor loading matrix
> ### Aliases: factor.rotate
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> #using the Harman 24 mental tests, rotate the 2nd and 3rd factors 45 degrees
> pc <- principal(Harman74.cor$cov,4,rotate=TRUE)
> pcr45 <- factor.rotate(pc,45,2,3)
> pcr90 <- factor.rotate(pcr45,45,2,3)
> print(factor.congruence(pc,pcr45),digits=3) #poor congruence with original
    [,1] [,2]  [,3] [,4]
PC1    1 0.00  0.00    0
PC2    0 0.74 -0.74    0
PC3    0 0.67  0.67    0
PC4    0 0.00  0.00    1
> print(factor.congruence(pc,pcr90),digits=3) #factor 2 and 3 have been exchanged and 3 flipped
    [,1] [,2] [,3] [,4]
PC1    1    0    0    0
PC2    0    0   -1    0
PC3    0    1    0    0
PC4    0    0    0    1
> 
> 
> 
> cleanEx(); nameEx("factor.stats")
> ### * factor.stats
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor.stats
> ### Title: Find various goodness of fit statistics for factor analysis and
> ###   principal components
> ### Aliases: factor.stats factor.scores
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> v9 <- sim.hierarchical()
Loading required package: MASS
> f3 <- factor.minres(v9,3)
> factor.stats(v9,f3,n.obs=500)
Call: factor.stats(r = v9, f = f3, n.obs = 500)

Test of the hypothesis that 3 factors are sufficient.

The degrees of freedom for the model is 12  and the fit was  0 
The number of observations was  500  with Chi Square =  0  with prob <  1 

Measures of factor score adequacy             
 Correlation of scores with factors            0.78 0.73 0.67
Multiple R square of scores with factors       0.61 0.53 0.45
Minimum correlation of factor score estimates  0.22 0.07 -0.1 
> f3o <- factor.pa(v9,3,rotate="Promax")
> factor.stats(v9,f3o,n.obs=500)
Call: factor.stats(r = v9, f = f3o, n.obs = 500)

Test of the hypothesis that 3 factors are sufficient.

The degrees of freedom for the model is 12  and the fit was  0.44 
The number of observations was  500  with Chi Square =  216.22  with prob <  1.4e-39 

Measures of factor score adequacy             
 Correlation of scores with factors            1 0.91 0.78
Multiple R square of scores with factors       0.99 0.83 0.61
Minimum correlation of factor score estimates  0.98 0.66 0.21 
> 
> 
> 
> 
> cleanEx(); nameEx("factor2cluster")
> ### * factor2cluster
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: factor2cluster
> ### Title: Extract cluster definitions from factor loadings
> ### Aliases: factor2cluster
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> ## Not run: 
> ##D f  <- factanal(x,4,covmat=Harman74.cor$cov)
> ##D factor2cluster(f) 
> ## End(Not run)
> #                       Factor1 Factor2 Factor3 Factor4
> #VisualPerception             0       1       0       0
> #Cubes                        0       1       0       0
> #PaperFormBoard               0       1       0       0
> #Flags                        0       1       0       0
> #GeneralInformation           1       0       0       0
> #PargraphComprehension        1       0       0       0
> #SentenceCompletion           1       0       0       0
> #WordClassification           1       0       0       0
> #WordMeaning                  1       0       0       0
> #Addition                     0       0       1       0
> #Code                         0       0       1       0
> #CountingDots                 0       0       1       0
> #StraightCurvedCapitals       0       0       1       0
> #WordRecognition              0       0       0       1
> #NumberRecognition            0       0       0       1
> #FigureRecognition            0       0       0       1
> #ObjectNumber                 0       0       0       1
> #NumberFigure                 0       0       0       1
> #FigureWord                   0       0       0       1
> #Deduction                    0       1       0       0
> #NumericalPuzzles             0       0       1       0
> #ProblemReasoning             0       1       0       0
> #SeriesCompletion             0       1       0       0
> #ArithmeticProblems           0       0       1       0
> 
> 
> 
> 
> 
> 
> cleanEx(); nameEx("fisherz")
> ### * fisherz
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fisherz
> ### Title: Fisher r to z and z to r and confidence intervals
> ### Aliases: fisherz fisherz2r r.con r2t
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> cors <- seq(-.9,.9,.1)
> zs <- fisherz(cors)
> rs <- fisherz2r(zs)
> round(zs,2)
 [1] -1.47 -1.10 -0.87 -0.69 -0.55 -0.42 -0.31 -0.20 -0.10  0.00  0.10  0.20
[13]  0.31  0.42  0.55  0.69  0.87  1.10  1.47
>  n <- 30
>  r <- seq(0,.9,.1)
>  rc <- matrix(r.con(r,n),ncol=2)
> t <- r*sqrt(n-2)/sqrt(1-r^2)
>  p <- (1-pt(t,n-2))/2
>  r.rc <- r.rc <- data.frame(r=r,z=fisherz(r),lower=rc[,1],upper=rc[,2],t=t,p=p)
>  round(r.rc,2)
     r    z lower upper     t    p
1  0.0 0.00 -0.36  0.36  0.00 0.25
2  0.1 0.10 -0.27  0.44  0.53 0.15
3  0.2 0.20 -0.17  0.52  1.08 0.07
4  0.3 0.31 -0.07  0.60  1.66 0.03
5  0.4 0.42  0.05  0.66  2.31 0.01
6  0.5 0.55  0.17  0.73  3.06 0.00
7  0.6 0.69  0.31  0.79  3.97 0.00
8  0.7 0.87  0.45  0.85  5.19 0.00
9  0.8 1.10  0.62  0.90  7.06 0.00
10 0.9 1.47  0.80  0.95 10.93 0.00
> 
> 
> 
> 
> cleanEx(); nameEx("galton")
> ### * galton
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: galton
> ### Title: Galton's Mid parent child height data
> ### Aliases: galton
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(galton)
> describe(galton)
       var   n  mean   sd median trimmed  mad  min  max range  skew kurtosis
parent   1 928 68.31 1.79   68.5   68.32 1.48 64.0 73.0     9 -0.04     0.05
child    2 928 68.09 2.52   68.2   68.12 2.97 61.7 73.7    12 -0.09    -0.35
         se
parent 0.06
child  0.08
> pairs.panels(galton,lm=TRUE)
> 
> 
> 
> cleanEx(); nameEx("geometric.mean")
> ### * geometric.mean
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: geometric.mean
> ### Title: Find the geometric mean of a vector or columns of a data.frame.
> ### Aliases: geometric.mean
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> 
> x <- seq(1,5)
> x2 <- x^2
> geometric.mean(x)
[1] 2.605171
> geometric.mean(x2)
[1] 6.786916
> 
> 
> 
> 
> cleanEx(); nameEx("guttman")
> ### * guttman
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: guttman
> ### Title: Alternative estimates of test reliabiity
> ### Aliases: tenberge glb guttman
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(attitude)
> glb(attitude)
$beta
[1] 0.7

$beta.factor
[1] 0.6

$alpha.pc
[1] 0.7309206

$glb.max
[1] 0.91

$glb.IC
[1] 0.88

$glb.Km
[1] 0.91

$glb.Fa
[1] 0.9

$r.smc
[1] 0.8752151

$tenberge
$tenberge$mu0
[1] 0.84

$tenberge$mu1
[1] 0.85

$tenberge$mu2
[1] 0.85

$tenberge$mu3
[1] 0.85


$keys
           IC1 IC2 ICr1 ICr2 K1 K2 F1 F2 f1 f2
rating       1   0    1    0  0  1  0  1  1  0
complaints   1   0    0    1  1  0  0  1  1  0
privileges   1   0    1    0  0  1  1  0  1  0
learning     1   0    0    1  0  1  1  0  1  0
raises       0   1    1    0  0  1  1  0  1  0
critical     0   1    1    0  0  1  1  0  0  1
advance      0   1    1    0  1  0  0  1  0  1

> guttman(attitude)
Call: guttman(r = attitude)

Alternative estimates of reliability
Beta =  0.7  This is an estimate of the worst split half reliability
Guttman bounds 
L1 =  0.72 
L2 =  0.85 
L3 (alpha) =  0.84 
L4 (max) =  0.91 
L5 =  0.83 
L6 (smc) =  0.88 
TenBerge bounds 
mu0 =  mu1 =  0.85 mu2 =  0.85 mu3 =  0.85 

alpha of first PC =  0.85 
estimated greatest lower bound =  0.91 

beta estimated by first and second PC =  0.62  This is an exploratory statistic 
> 
> 
> 
> 
> cleanEx(); nameEx("harmonic.mean")
> ### * harmonic.mean
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: harmonic.mean
> ### Title: Find the harmonic mean of a vector, matrix, or columns of a
> ###   data.frame
> ### Aliases: harmonic.mean
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> x <- seq(1,5)
> x2 <- x^2
> harmonic.mean(x)
[1] 2.189781
> harmonic.mean(x2)
[1] 3.416208
> 
> 
> 
> 
> cleanEx(); nameEx("headtail")
> ### * headtail
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: headtail
> ### Title: Combine calls to head and tail
> ### Aliases: headtail
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> x <- matrix(sample(10,1000,TRUE),ncol=5)
> headtail(x,4,8)
     X1  X2  X3  X4  X5
1     3   3   7   9   9
2     4   3   2  10   1
3     6   6  10   2  10
4    10   3   9   8   8
... ... ... ... ... ...
193   2  10   8  10   7
194  10   6   7   4   3
195   3   2   9   6   2
196   6   4   7   7   8
197   2   7   5   3   7
198   9   4   5   6   3
199   4   4   4   9   2
200   8   2   5   4   3
> 
> 
> 
> 
> cleanEx(); nameEx("heights")
> ### * heights
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: heights
> ### Title: A data.frame of the Galton (1888) height and cubit data set.
> ### Aliases: heights
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(heights)
> ellipses(heights,n=1,main="Galton's co-relation data set")
> 
> 
> 
> 
> cleanEx(); nameEx("interp.median")
> ### * interp.median
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: interp.median
> ### Title: Find the interpolated sample median, quartiles, or specific
> ###   quantiles for a vector, matrix, or data frame
> ### Aliases: interp.median interp.quantiles interp.quartiles interp.boxplot
> ###   interp.values interp.qplot.by interp.q interp.quart
> ### Keywords: univar
> 
> ### ** Examples
> 
> interp.median(c(1,2,3,3,3))  # compare with median = 3
[1] 2.666667
> interp.median(c(1,2,2,5))
[1] 2
> interp.quantiles(c(1,2,2,5),.25)
> x <- sample(10,100,TRUE)
> interp.quartiles(x)
      Q1   Median       Q3 
3.571429 5.357143 7.966667 
> #
> x <-  c(1,1,2,2,2,3,3,3,3,4,5,1,1,1,2,2,3,3,3,3,4,5,1,1,1,2,2,3,3,3,3,4,2)
> y <-  c(1,2,3,3,3,3,4,4,4,5,5,1,2,3,3,3,3,4,4,5,5,5,1,5,3,3,3,3,4,4,4,5,5)
> x <-  x[order(x)]   #sort the data by ascending order to make it clearer
> y <- y[order(y)]
> xv <- interp.values(x)
> yv <- interp.values(y)
> barplot(x,space=0,xlab="ordinal position",ylab="value")
> lines(1:length(x)-.5,xv)
> points(c(length(x)/4,length(x)/2,3*length(x)/4),interp.quartiles(x))
> barplot(y,space=0,xlab="ordinal position",ylab="value")
> lines(1:length(y)-.5,yv)
> points(c(length(y)/4,length(y)/2,3*length(y)/4),interp.quartiles(y))
> data(galton)
> interp.median(galton)
       parent    child
[1,] 68.32877 68.18333
> interp.qplot.by(galton$child,galton$parent,ylab="child height"
+ ,xlab="Mid parent height") 
> 
> 
> 
> 
> cleanEx(); nameEx("iqitems")
> ### * iqitems
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: iqitems
> ### Title: 14 multiple choice IQ items
> ### Aliases: iqitems
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(iqitems)
> iq.keys <- c(4,4,3,1,4,3,2,3,1,4,1,3,4,3)
> score.multiple.choice(iq.keys,iqitems)
$item.stats
     key    0    1    2    3    4    5    6    r    n mean   sd  skew kurtosis
iq1    4 0.04 0.01 0.03 0.09 0.80 0.02 0.01 0.59 1000 0.80 0.40 -1.51     0.27
iq8    4 0.03 0.10 0.01 0.02 0.80 0.01 0.04 0.39 1000 0.80 0.40 -1.49     0.22
iq10   3 0.10 0.22 0.09 0.37 0.04 0.13 0.04 0.35 1000 0.37 0.48  0.53    -1.72
iq15   1 0.03 0.65 0.16 0.15 0.00 0.00 0.00 0.35 1000 0.65 0.48 -0.63    -1.60
iq20   4 0.03 0.02 0.03 0.03 0.85 0.02 0.01 0.42 1000 0.85 0.35 -2.00     2.01
iq44   3 0.03 0.10 0.06 0.64 0.02 0.14 0.01 0.42 1000 0.64 0.48 -0.61    -1.64
iq47   2 0.04 0.08 0.59 0.06 0.11 0.07 0.05 0.51 1000 0.59 0.49 -0.35    -1.88
iq2    3 0.07 0.08 0.31 0.32 0.15 0.05 0.02 0.26 1000 0.32 0.46  0.80    -1.37
iq11   1 0.04 0.87 0.03 0.01 0.01 0.01 0.04 0.54 1000 0.87 0.34 -2.15     2.61
iq16   4 0.05 0.05 0.08 0.07 0.74 0.01 0.00 0.56 1000 0.74 0.44 -1.11    -0.77
iq32   1 0.04 0.54 0.02 0.14 0.10 0.04 0.12 0.50 1000 0.54 0.50 -0.17    -1.97
iq37   3 0.07 0.10 0.09 0.26 0.13 0.02 0.34 0.23 1000 0.26 0.44  1.12    -0.74
iq43   4 0.04 0.07 0.04 0.02 0.78 0.03 0.00 0.50 1000 0.78 0.41 -1.35    -0.18
iq49   3 0.06 0.27 0.09 0.32 0.14 0.08 0.05 0.28 1000 0.32 0.47  0.79    -1.38
       se
iq1  0.01
iq8  0.01
iq10 0.02
iq15 0.02
iq20 0.01
iq44 0.02
iq47 0.02
iq2  0.01
iq11 0.01
iq16 0.01
iq32 0.02
iq37 0.01
iq43 0.01
iq49 0.01

$alpha
         Averages
Averages     0.63

$av.r
         Averages
Averages     0.11

> 
> 
> 
> cleanEx(); nameEx("kappa")
> ### * kappa
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: wkappa
> ### Title: Find Cohen's kappa and weighted kappa coefficients for
> ###   correlation of two raters
> ### Aliases: wkappa
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> 
> cohen <- matrix(c(
+ 0.44, 0.05, 0.01,
+ 0.07, 0.20, 0.03,
+ 0.09, 0.05, 0.06),ncol=3)
> 
> wkappa(cohen)
$kappa
[1] 0.4915254

$weighted.kappa
[1] 0.4814815

> 
> fleiss <- matrix(c(
+  0.53, 0.05, 0.02,
+  0.11, 0.14, 0.05,
+  0.01, 0.06, 0.03),ncol=3)
> 
> weights <- matrix(c(
+  1.0000, 0.0000, 0.4444,
+  0.0000, 1.0000, 0.6666,
+  0.4444, 0.6666, 1.0000),ncol=3)
>  
>  wkappa(fleiss,weights)
$kappa
[1] 0.4285714

$weighted.kappa
[1] 0.5070508

> 
> 
> 
> cleanEx(); nameEx("make.keys")
> ### * make.keys
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: make.keys
> ### Title: Create a keys matrix for use by score.items or cluster.cor
> ### Aliases: make.keys
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> data(attitude)
>  key.list <- list(all=c(1,2,3,4,-5,6,7),
+                   first=c(1,2,3),
+                   last=c(4,5,6,7))
>  keys <- make.keys(7,key.list,item.labels = colnames(attitude))
>  keys
           all first last
rating       1     1    0
complaints   1     1    0
privileges   1     1    0
learning     1     0    1
raises      -1     0    1
critical     1     0    1
advance      1     0    1
>  
>  scores <- score.items(keys,attitude,short=TRUE)
>  scores
Call: score.items(keys = keys, items = attitude, short = TRUE)

(Unstandardized) Alpha:
       all first last
alpha 0.54  0.82 0.75

Average item correlation:
           all first last
average.r 0.14  0.61 0.42

 Guttman 6* reliability: 
          all first last
Lambda.6 0.77  0.84 0.79

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
       all first last
all   0.54  1.37 1.28
first 0.92  0.82 0.78
last  0.81  0.61 0.75

Item by scale correlations:
 corrected for item overlap and scale reliability
            all first last
rating     0.81  0.85 0.59
complaints 0.86  0.92 0.64
privileges 0.64  0.58 0.54
learning   0.79  0.73 0.73
raises     0.87  0.73 0.85
critical   0.23  0.21 0.36
advance    0.48  0.31 0.72
>  
>  data(bfi)
>  keys.list <- list(agree=c(-1,2:5),conscientious=c(6:8,-9,-10),extraversion=c(-11,-12,13:15),neuroticism=c(16:20),openness = c(21,-22,23,24,-25))
>  keys <- make.keys(25,keys.list,item.labels=colnames(bfi))
>  scores <- score.items(keys,bfi,short=TRUE)
>  scores
Call: score.items(keys = keys, items = bfi, short = TRUE)

(Unstandardized) Alpha:
      agree conscientious extraversion neuroticism openness
alpha  0.65          0.74         0.77         0.8     0.47

Average item correlation:
          agree conscientious extraversion neuroticism openness
average.r  0.27          0.36         0.41        0.46     0.15

 Guttman 6* reliability: 
         agree conscientious extraversion neuroticism openness
Lambda.6  0.66          0.74         0.78        0.82     0.53

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
              agree conscientious extraversion neuroticism openness
agree          0.65          0.36         0.62       -0.30     0.33
conscientious  0.25          0.74         0.28       -0.22     0.24
extraversion   0.44          0.21         0.77       -0.24     0.57
neuroticism   -0.22         -0.17        -0.19        0.81     0.11
openness       0.18          0.14         0.34        0.07     0.47

Item by scale correlations:
 corrected for item overlap and scale reliability
   agree conscientious extraversion neuroticism openness
A1 -0.40         -0.12        -0.07        0.16    -0.16
A2  0.63          0.20         0.41       -0.09     0.26
A3  0.63          0.20         0.41       -0.13     0.30
A4  0.39          0.21         0.22       -0.17    -0.05
A5  0.61          0.19         0.54       -0.23     0.17
C1  0.12          0.59         0.11       -0.03     0.23
C2  0.19          0.56         0.08        0.07     0.14
C3  0.24          0.61         0.23       -0.03     0.15
C4 -0.29         -0.70        -0.22        0.30    -0.17
C5 -0.22         -0.57        -0.20        0.32    -0.01
E1 -0.36         -0.07        -0.65        0.09    -0.42
E2 -0.35         -0.21        -0.69        0.28    -0.23
E3  0.44          0.19         0.61       -0.09     0.50
E4  0.53          0.12         0.68       -0.23     0.10
E5  0.27          0.34         0.58       -0.05     0.49
N1 -0.26         -0.17        -0.06        0.77     0.13
N2 -0.26         -0.11        -0.06        0.74     0.09
N3 -0.19         -0.16        -0.10        0.71     0.10
N4 -0.24         -0.27        -0.37        0.63     0.08
N5 -0.06         -0.04        -0.23        0.54    -0.03
O1  0.21          0.17         0.30       -0.01     0.54
O2  0.13          0.13        -0.22       -0.06    -0.09
O3  0.29          0.21         0.43        0.01     0.63
O4  0.17          0.10         0.00        0.22     0.42
O5 -0.12         -0.14        -0.15        0.06    -0.42
> 
> 
> 
> 
> cleanEx(); nameEx("mat.regress")
> ### * mat.regress
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mat.regress
> ### Title: Multiple Regression from matrix input
> ### Aliases: mat.regress
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> 
> ## Not run: 
> ##D test.data <- Harman74.cor$cov     #24 mental variables
> ##D #choose 3 of them to regress against another 4 -- arbitrary choice of variables
> ##D print(mat.regress(test.data,c(1,2,3),c(4,5,10,12)),digits=2)
> ## End(Not run)
> #gives this output
> #print(mat.regress(test.data,c(1,2,3),c(4,5,10,12)),digits=2)
> #$beta
> #                 Flags GeneralInformation Addition CountingDots
> #VisualPerception  0.40               0.22     0.16         0.30
> #Cubes             0.06               0.18     0.06         0.05
> #PaperFormBoard    0.12               0.10    -0.16         0.00
> #
> #$R
> #             Flags GeneralInformation           Addition       CountingDots 
> #             0.49               0.38               0.18               0.32 
> #
> #$R2
> #             Flags GeneralInformation           Addition       CountingDots 
> #              0.24               0.15               0.03               0.10 
> #
> #
> data(attitude)
> mat.regress(attitude,c(1:3),c(4:7)) #standardized regression from raw data
$beta
           learning raises critical advance
rating         0.42   0.11     0.01   -0.05
complaints     0.08   0.39     0.11    0.08
privileges     0.24   0.09     0.05    0.26

$se
           learning raises critical advance
rating         0.25   0.22     0.28    0.28
complaints     0.25   0.22     0.28    0.28
privileges     0.17   0.15     0.19    0.19

$t
           learning raises critical advance
rating         1.71   0.51     0.03   -0.19
complaints     0.33   1.79     0.39    0.28
privileges     1.46   0.64     0.27    1.41

$Probability
           learning raises critical advance
rating         0.10   0.61     0.97    0.85
complaints     0.75   0.09     0.70    0.78
privileges     0.16   0.53     0.79    0.17

$R
learning   raises critical  advance 
    0.67     0.68     0.19     0.35 

$R2
learning   raises critical  advance 
    0.45     0.46     0.04     0.12 

$shrunkenR2
learning   raises critical  advance 
    0.39     0.40    -0.07     0.02 

$seR2
learning   raises critical  advance 
    0.11     0.11     0.06     0.09 

$F
learning   raises critical  advance 
    7.22     7.40     0.34     1.19 

$probF
learning   raises critical  advance 
   0.001    0.001    0.796    0.333 

$df
[1]  3 26

> 
> 
> 
> cleanEx(); nameEx("mat.sort")
> ### * mat.sort
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mat.sort
> ### Title: Sort the elements of a correlation matrix to reflect factor
> ###   loadings
> ### Aliases: mat.sort
> ### Keywords: multivariate models
> 
> ### ** Examples
> data(bifactor)
> sorted <- mat.sort(Bechtoldt.1,fa(Bechtoldt.1,5))
> cor.plot(sorted)
> 
> 
> 
> cleanEx(); nameEx("matrix.addition")
> ### * matrix.addition
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: matrix.addition
> ### Title: A function to add two vectors or matrices
> ### Aliases: matrix.addition %+%
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> 
> x <- seq(1,4)
> z <- x %+% -t(x)
> x
[1] 1 2 3 4
> z
     [,1] [,2] [,3] [,4]
[1,]    0   -1   -2   -3
[2,]    1    0   -1   -2
[3,]    2    1    0   -1
[4,]    3    2    1    0
> #compare with outer(x,-x,FUN="+")
> x <- matrix(seq(1,6),ncol=2)
> y <- matrix(seq(1,10),nrow=2)
> z <- x %+% y
> x
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
> y
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    3    5    7    9
[2,]    2    4    6    8   10
> z
     [,1] [,2] [,3] [,4] [,5]
[1,]    8   12   16   20   24
[2,]   10   14   18   22   26
[3,]   12   16   20   24   28
> #but compare this with outer(x ,y,FUN="+") 
> 
> 
> 
> cleanEx(); nameEx("msq")
> ### * msq
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: msq
> ### Title: 75 mood items from the Motivational State Questionnaire for 3896
> ###   participants
> ### Aliases: msq
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(msq)
> describe(msq)
             var    n   mean     sd median trimmed   mad  min     max  range
MSQ_Time       1 3895  13.49   4.27  13.23   13.22  5.78 5.75   22.75   17.0
active         2 3890   1.03   0.93   1.00    0.95  1.48 0.00    3.00    3.0
afraid         3 3891   0.12   0.40   0.00    0.00  0.00 0.00    3.00    3.0
alert          4 3885   1.15   0.91   1.00    1.09  1.48 0.00    3.00    3.0
alone          5 1838   0.54   0.77   0.00    0.38  0.00 0.00    3.00    3.0
angry          6 3887   0.22   0.56   0.00    0.07  0.00 0.00    3.00    3.0
aroused        7 3890   0.71   0.85   0.00    0.59  0.00 0.00    3.00    3.0
ashamed        8 3885   0.12   0.43   0.00    0.00  0.00 0.00    3.00    3.0
astonished     9 3883   0.13   0.41   0.00    0.00  0.00 0.00    3.00    3.0
at-ease       10 3879   1.59   0.92   2.00    1.61  1.48 0.00    3.00    3.0
at-rest       11 3879   1.20   0.92   1.00    1.13  1.48 0.00    3.00    3.0
attentive     12 3890   1.18   0.86   1.00    1.14  1.48 0.00    3.00    3.0
blue          13 3891   0.37   0.65   0.00    0.24  0.00 0.00    3.00    3.0
bored         14 3892   1.31   0.94   1.00    1.26  1.48 0.00    3.00    3.0
calm          15 3814   1.55   0.92   2.00    1.56  1.48 0.00    3.00    3.0
clutched-up   16 3873   0.37   0.66   0.00    0.23  0.00 0.00    3.00    3.0
confident     17 3889   1.49   0.92   2.00    1.49  1.48 0.00    3.00    3.0
content       18 3874   1.45   0.93   1.00    1.43  1.48 0.00    3.00    3.0
delighted     19 3890   0.78   0.84   1.00    0.68  1.48 0.00    3.00    3.0
depressed     20 3879   0.41   0.68   0.00    0.27  0.00 0.00    3.00    3.0
determined    21 3889   1.04   0.95   1.00    0.95  1.48 0.00    3.00    3.0
distressed    22 3888   0.40   0.71   0.00    0.23  0.00 0.00    3.00    3.0
drowsy        23 3884   1.16   1.03   1.00    1.08  1.48 0.00    3.00    3.0
dull          24 3887   0.76   0.85   1.00    0.64  1.48 0.00    3.00    3.0
elated        25 3881   0.52   0.76   0.00    0.38  0.00 0.00    3.00    3.0
energetic     26 3890   0.84   0.89   1.00    0.74  1.48 0.00    3.00    3.0
enthusiastic  27 3890   0.92   0.90   1.00    0.83  1.48 0.00    3.00    3.0
excited       28 3890   0.72   0.86   0.00    0.60  0.00 0.00    3.00    3.0
fearful       29 3876   0.13   0.43   0.00    0.00  0.00 0.00    3.00    3.0
frustrated    30 3885   0.49   0.80   0.00    0.32  0.00 0.00    3.00    3.0
full-of-pep   31 3884   0.81   0.91   1.00    0.69  1.48 0.00    3.00    3.0
gloomy        32 3884   0.40   0.67   0.00    0.26  0.00 0.00    3.00    3.0
grouchy       33 3891   0.47   0.73   0.00    0.31  0.00 0.00    3.00    3.0
guilty        34 3891   0.18   0.50   0.00    0.05  0.00 0.00    3.00    3.0
happy         35 3880   1.22   0.95   1.00    1.16  1.48 0.00    3.00    3.0
hostile       36 3885   0.29   0.60   0.00    0.16  0.00 0.00    3.00    3.0
inspired      37 3890   0.59   0.82   0.00    0.44  0.00 0.00    3.00    3.0
intense       38 3889   0.69   0.86   0.00    0.56  0.00 0.00    3.00    3.0
interested    39 3884   1.14   0.90   1.00    1.09  1.48 0.00    3.00    3.0
irritable     40 3880   0.55   0.77   0.00    0.41  0.00 0.00    3.00    3.0
jittery       41 3890   0.59   0.80   0.00    0.45  0.00 0.00    3.00    3.0
kindly        42 1836   1.18   0.91   1.00    1.14  1.48 0.00    3.00    3.0
lively        43 3886   0.92   0.89   1.00    0.84  1.48 0.00    3.00    3.0
lonely        44 3890   0.47   0.73   0.00    0.31  0.00 0.00    3.00    3.0
nervous       45 3879   0.35   0.65   0.00    0.22  0.00 0.00    3.00    3.0
placid        46 3877   1.20   0.90   1.00    1.15  1.48 0.00    3.00    3.0
pleased       47 3883   0.90   0.87   1.00    0.82  1.48 0.00    3.00    3.0
proud         48 3889   0.94   0.95   1.00    0.84  1.48 0.00    3.00    3.0
quiescent     49 3760   0.65   0.73   1.00    0.55  1.48 0.00    3.00    3.0
quiet         50 3891   1.39   0.94   1.00    1.37  1.48 0.00    3.00    3.0
relaxed       51 3889   1.68   0.88   2.00    1.72  1.48 0.00    3.00    3.0
sad           52 3886   0.34   0.63   0.00    0.22  0.00 0.00    3.00    3.0
satisfied     53 3889   1.34   0.89   1.00    1.32  1.48 0.00    3.00    3.0
scared        54 3886   0.17   0.48   0.00    0.04  0.00 0.00    3.00    3.0
scornful      55 1838   0.17   0.46   0.00    0.05  0.00 0.00    3.00    3.0
serene        56 3884   1.18   0.89   1.00    1.12  1.48 0.00    3.00    3.0
sleepy        57 3880   1.25   1.05   1.00    1.18  1.48 0.00    3.00    3.0
sluggish      58 3888   1.17   0.97   1.00    1.09  1.48 0.00    3.00    3.0
sociable      59 3890   1.34   0.89   1.00    1.31  1.48 0.00    3.00    3.0
sorry         60 3881   0.20   0.55   0.00    0.06  0.00 0.00    3.00    3.0
still         61 3884   1.18   0.89   1.00    1.12  1.48 0.00    3.00    3.0
strong        62 3889   1.05   0.93   1.00    0.98  1.48 0.00    3.00    3.0
surprised     63 3890   0.16   0.45   0.00    0.03  0.00 0.00    3.00    3.0
tense         64 3886   0.52   0.77   0.00    0.37  0.00 0.00    3.00    3.0
tired         65 3886   1.39   1.04   1.00    1.36  1.48 0.00    3.00    3.0
unhappy       66 3891   0.37   0.65   0.00    0.24  0.00 0.00    3.00    3.0
upset         67 3888   0.31   0.63   0.00    0.17  0.00 0.00    3.00    3.0
vigorous      68 3886   0.64   0.81   0.00    0.52  0.00 0.00    3.00    3.0
wakeful       69 3886   1.19   0.90   1.00    1.14  1.48 0.00    3.00    3.0
warmhearted   70 3889   1.45   0.93   1.00    1.44  1.48 0.00    3.00    3.0
wide-awake    71 3884   0.94   0.95   1.00    0.83  1.48 0.00    3.00    3.0
anxious       72 2047   0.67   0.86   0.00    0.54  0.00 0.00    3.00    3.0
idle          73 2050   0.96   0.95   1.00    0.85  1.48 0.00    9.00    9.0
cheerful      74 2046   0.98   0.91   1.00    0.90  1.48 0.00    3.00    3.0
inactive      75 2050   1.08   0.97   1.00    0.97  1.48 0.00    3.00    3.0
tranquil      76 2053   1.29   0.87   1.00    1.24  1.48 0.00    3.00    3.0
EA            77 3298  11.74   7.43  11.00   11.42  7.41 0.00   30.00   30.0
TA            78 3255  11.34   4.95  11.00   11.13  4.45 0.00   33.00   33.0
PA            79 3324   9.73   6.75   9.00    9.29  7.41 0.00   30.00   30.0
NegAff        80 3304   2.98   3.80   2.00    2.25  2.97 0.00   28.00   28.0
Extraversion  81 3890  13.55   4.12  14.00   13.75  4.45 0.00   24.00   24.0
Neuroticism   82 3890  10.39   4.90  10.00   10.27  5.93 0.00   24.00   24.0
Lie           83 3890   2.11   1.52   2.00    2.01  1.48 0.00    9.00    9.0
Sociability   84 3890   7.72   2.81   8.00    7.94  2.97 0.00   13.00   13.0
Impulsivity   85 3890   4.45   1.87   4.00    4.43  1.48 0.00    9.00    9.0
MSQ_Round     86 2797  13.44   4.39  13.00   13.15  5.93 6.00   23.00   17.0
scale*        87 3896   1.53   0.50   2.00    1.53  0.00 1.00    2.00    1.0
ID            88 3896 159.38 363.77  67.00   85.40 68.20 1.00 4184.00 4183.0
exper*        89 3896  21.10  10.48  23.00   21.31 11.86 1.00   39.00   38.0
condition     90 3896   2.57   1.50   2.00    2.46  1.48 1.00    5.00    4.0
TOD           91 3895  13.19   4.26  13.00   12.93  5.93 5.00   22.50   17.5
TOD24         92 2794  18.13  10.76  15.00   16.45  8.15 5.00   46.50   41.5
              skew kurtosis   se
MSQ_Time      0.38    -1.19 0.07
active        0.47    -0.76 0.01
afraid        4.18    19.99 0.01
alert         0.33    -0.76 0.01
alone         1.41     1.43 0.02
angry         2.92     8.85 0.01
aroused       0.95    -0.04 0.01
ashamed       4.23    19.94 0.01
astonished    3.62    14.21 0.01
at-ease      -0.09    -0.83 0.01
at-rest       0.33    -0.74 0.01
attentive     0.29    -0.61 0.01
blue          1.79     2.94 0.01
bored         0.30    -0.78 0.02
calm         -0.01    -0.83 0.01
clutched-up   1.79     2.72 0.01
confident    -0.05    -0.85 0.01
content       0.02    -0.86 0.01
delighted     0.72    -0.46 0.01
depressed     1.69     2.44 0.01
determined    0.48    -0.81 0.02
distressed    1.87     2.98 0.01
drowsy        0.46    -0.93 0.02
dull          0.91     0.04 0.01
elated        1.33     1.00 0.01
energetic     0.74    -0.42 0.01
enthusiastic  0.60    -0.58 0.01
excited       0.96    -0.01 0.01
fearful       3.81    16.05 0.01
frustrated    1.57     1.60 0.01
full-of-pep   0.83    -0.34 0.01
gloomy        1.74     2.71 0.01
grouchy       1.52     1.71 0.01
guilty        3.14    10.48 0.01
happy         0.23    -0.91 0.02
hostile       2.15     4.39 0.01
inspired      1.25     0.73 0.01
intense       1.03     0.11 0.01
interested    0.32    -0.73 0.01
irritable     1.30     1.02 0.01
jittery       1.24     0.81 0.01
kindly        0.20    -0.89 0.02
lively        0.59    -0.58 0.01
lonely        1.54     1.84 0.01
nervous       1.93     3.47 0.01
placid        0.27    -0.76 0.01
pleased       0.57    -0.59 0.01
proud         0.60    -0.75 0.02
quiescent     0.86     0.17 0.01
quiet         0.15    -0.85 0.02
relaxed      -0.17    -0.68 0.01
sad           1.96     3.79 0.01
satisfied     0.00    -0.83 0.01
scared        3.25    11.36 0.01
scornful      3.16    11.20 0.01
serene        0.35    -0.64 0.01
sleepy        0.40    -1.04 0.02
sluggish      0.46    -0.74 0.02
sociable      0.03    -0.81 0.01
sorry         3.07     9.76 0.01
still         0.33    -0.66 0.01
strong        0.40    -0.87 0.01
surprised     3.24    11.16 0.01
tense         1.41     1.34 0.01
tired         0.22    -1.10 0.02
unhappy       1.89     3.48 0.01
upset         2.27     5.14 0.01
vigorous      1.03     0.09 0.01
wakeful       0.27    -0.74 0.01
warmhearted   0.00    -0.88 0.01
wide-awake    0.65    -0.63 0.02
anxious       1.09     0.26 0.02
idle          1.17     4.26 0.02
cheerful      0.51    -0.73 0.02
inactive      0.49    -0.80 0.02
tranquil      0.24    -0.61 0.02
EA            0.34    -0.66 0.13
TA            0.48     0.27 0.09
PA            0.51    -0.42 0.12
NegAff        2.30     7.31 0.07
Extraversion -0.41    -0.22 0.07
Neuroticism   0.19    -0.59 0.08
Lie           0.68     0.31 0.02
Sociability  -0.62    -0.20 0.05
Impulsivity   0.07    -0.56 0.03
MSQ_Round     0.41    -1.10 0.08
scale*       -0.11    -1.99 0.01
ID            5.75    40.62 5.83
exper*       -0.21    -1.02 0.17
condition     0.65    -1.10 0.02
TOD           0.37    -1.21 0.07
TOD24         1.17     0.11 0.20
> 
> 
> 
> cleanEx(); nameEx("multi.hist")
> ### * multi.hist
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: multi.hist
> ### Title: Multiple histograms with density and normal fits on one page
> ### Aliases: multi.hist histo.density
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> #multi.hist(attitude[-1]) 
> 
> 
> 
> cleanEx(); nameEx("neo")
> ### * neo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: neo
> ### Title: NEO correlation matrix from the NEO_PI_R manual
> ### Aliases: neo
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(neo)
> n5 <- factor.minres(neo,5)
> neo.keys <- make.keys(30,list(N=c(1:6),E=c(7:12),O=c(13:18),A=c(19:24),C=c(25:30)))
> n5p <- target.rot(n5,neo.keys) #show a targeted rotation for simple structure
> n5p

Call: NULL
    V   MR1   MR4   MR3   MR2   MR5   h2   u2
N1  1  0.81                         0.67 0.33
N2  2  0.62             -0.40       0.55 0.45
N3  3  0.79                         0.65 0.35
N4  4  0.69                         0.52 0.48
N5  5  0.43  0.34                   0.38 0.62
N6  6  0.63                         0.48 0.52
E1  7        0.59        0.46       0.57 0.43
E2  8        0.57                   0.38 0.62
E3  9        0.45                   0.37 0.63
E4 10        0.53              0.35 0.46 0.54
E5 11        0.56                   0.39 0.61
E6 12        0.71                   0.56 0.44
O1 13              0.50             0.36 0.64
O2 14              0.70             0.56 0.44
O3 15  0.39  0.33  0.42             0.45 0.55
O4 16              0.43             0.24 0.76
O5 17              0.71             0.57 0.43
O6 18              0.33             0.15 0.85
A1 19 -0.30              0.49       0.37 0.63
A2 20                    0.57       0.44 0.56
A3 21        0.48        0.61       0.64 0.36
A4 22                    0.69       0.53 0.47
A5 23                    0.50       0.32 0.68
A6 24                    0.58       0.38 0.62
C1 25                          0.54 0.39 0.61
C2 26                          0.65 0.46 0.54
C3 27                          0.64 0.46 0.54
C4 28                          0.71 0.58 0.42
C5 29                          0.70 0.55 0.45
C6 30                          0.53 0.39 0.61

                MR1  MR4  MR3  MR2  MR5
SS loadings    3.32 2.90 1.88 2.78 2.96
Proportion Var 0.11 0.10 0.06 0.09 0.10
Cumulative Var 0.11 0.21 0.27 0.36 0.46

 With factor correlations of 
      MR1   MR4   MR3   MR2   MR5
MR1  1.00 -0.17 -0.05 -0.18 -0.47
MR4 -0.17  1.00  0.41 -0.04  0.24
MR3 -0.05  0.41  1.00 -0.03  0.04
MR2 -0.18 -0.04 -0.03  1.00  0.13
MR5 -0.47  0.24  0.04  0.13  1.00
> 
> 
> 
> 
> cleanEx(); nameEx("omega")
> ### * omega
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: omega
> ### Title: Calculate McDonald's omega estimates of general and total factor
> ###   saturation
> ### Aliases: omega
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D test.data <- Harman74.cor$cov
> ##D my.omega <- omega(test.data,3)       
> ##D print(my.omega,digits=2)
> ##D #
> ##D          
> ## End(Not run)
> #create 9 variables with a hierarchical structure
> jen.data <- sim.hierarchical()
Loading required package: MASS
> #with correlations of
> jen.data
       V1     V2     V3     V4     V5    V6     V7     V8     V9
V1 1.0000 0.5600 0.4800 0.4032 0.3456 0.288 0.3024 0.2520 0.2016
V2 0.5600 1.0000 0.4200 0.3528 0.3024 0.252 0.2646 0.2205 0.1764
V3 0.4800 0.4200 1.0000 0.3024 0.2592 0.216 0.2268 0.1890 0.1512
V4 0.4032 0.3528 0.3024 1.0000 0.4200 0.350 0.2352 0.1960 0.1568
V5 0.3456 0.3024 0.2592 0.4200 1.0000 0.300 0.2016 0.1680 0.1344
V6 0.2880 0.2520 0.2160 0.3500 0.3000 1.000 0.1680 0.1400 0.1120
V7 0.3024 0.2646 0.2268 0.2352 0.2016 0.168 1.0000 0.3000 0.2400
V8 0.2520 0.2205 0.1890 0.1960 0.1680 0.140 0.3000 1.0000 0.2000
V9 0.2016 0.1764 0.1512 0.1568 0.1344 0.112 0.2400 0.2000 1.0000
> #find omega 
> if(require(Rgraphviz)) {jen.omega <- omega(jen.data,digits=2)} else {jen.omega <- omega(jen.data,digits=2,plot=FALSE)}
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
Loading required package: GPArotation
> jen.omega
Omega 
Call: omega(m = jen.data, digits = 2)
Alpha:                 0.76 
G.6:                   0.76 
Omega Hierarchical:    0.69 
Omega H asymptotic:    0.86 
Omega Total            0.8 

Schmid Leiman Factor loadings greater than  0.2 
      g  F1*  F2*  F3*   h2   u2
V1 0.72 0.35           0.64 0.36
V2 0.63 0.31           0.49 0.51
V3 0.54 0.26           0.36 0.64
V4 0.56      0.42      0.49 0.51
V5 0.48      0.36      0.36 0.64
V6 0.40      0.30      0.25 0.75
V7 0.42           0.43 0.36 0.64
V8 0.35           0.36 0.25 0.75
V9 0.28           0.29      0.84

With eigenvalues of:
   g  F1*  F2*  F3* 
2.29 0.28 0.40 0.39 

general/max  5.78   max/min =   1.4
The degrees of freedom for the model is 12  and the fit was  0 

Measures of factor score adequacy              g F1* F2* F3*
 Correlation of scores with factors            0.85 0.46 0.57 0.57
Multiple R square of scores with factors       0.73 0.21 0.32 0.32
Minimum correlation of factor score estimates  0.46 -0.57 -0.35 -0.35 
> 
> #create 8 items with a two factor solution, showing the use of the flip option
> sim2 <- item.sim(8)
> omega(sim2)   #an example of misidentification-- remember to look at the loadings matrices.
Omega 
Call: omega(m = sim2)
Alpha:                 0.37 
G.6:                   0.51 
Omega Hierarchical:    0.56 
Omega H asymptotic:    0.87 
Omega Total            0.65 

Schmid Leiman Factor loadings greater than  0.2 
       g  F1*  F2*  F3*   h2   u2
V1       0.61           0.37 0.63
V2-      0.63           0.40 0.60
V3  0.53           0.20 0.32 0.68
V4  0.74                0.55 0.45
V5-      0.59           0.35 0.65
V6       0.65           0.44 0.56
V7- 0.49           0.87 1.00     
V8- 0.63                0.41 0.59

With eigenvalues of:
   g  F1*  F2*  F3* 
1.48 1.53 0.01 0.82 

general/max  0.97   max/min =   273.3
The degrees of freedom for the model is 7  and the fit was  0.01 
The number of observations was  500  with Chi Square =  5.24  with prob <  0.63 

Measures of factor score adequacy              g F1* F2* F3*
 Correlation of scores with factors            0.84 1.55 0.06 0.96
Multiple R square of scores with factors       0.71 2.41 0 0.93
Minimum correlation of factor score estimates  0.42 3.82 -0.99 0.86 
> omega(sim2,2)  #this shows that in fact there is no general factor
Three factors are required for identification -- general factor loadings set to be equal. Proceed with caution.
Three factors are required for identification -- general factor loadings set to be equal. Proceed with caution.
Omega 
Call: omega(m = sim2, nfactors = 2)
Alpha:                 0.63 
G.6:                   0.67 
Omega Hierarchical:    0.05 
Omega H asymptotic:    0.07 
Omega Total            0.72 

Schmid Leiman Factor loadings greater than  0.2 
       g  F1*  F2*   h2   u2
V1-           0.60 0.37 0.63
V2-           0.62 0.40 0.60
V3       0.60      0.37 0.63
V4       0.63      0.41 0.59
V5            0.57 0.34 0.66
V6            0.63 0.42 0.58
V7-      0.59      0.36 0.64
V8-      0.64      0.43 0.57

With eigenvalues of:
   g  F1*  F2* 
0.11 1.53 1.47 

general/max  0.07   max/min =   1.04
The degrees of freedom for the model is 13  and the fit was  0.06 
The number of observations was  500  with Chi Square =  31.63  with prob <  0.0027 

Measures of factor score adequacy              g F1* F2*
 Correlation of scores with factors            0.23 0.84 0.83
Multiple R square of scores with factors       0.05 0.7 0.69
Minimum correlation of factor score estimates  -0.9 0.4 0.38 
> omega(sim2,2,option="first") #but, if we define one of the two group factors as a general factor, we get a falsely high omega 
Three factors are required for identification -- general factor loading set to be 1 for group factor 1. Proceed with caution.
Three factors are required for identification -- general factor loading set to be 1 for group factor 1. Proceed with caution.
Omega 
Call: omega(m = sim2, nfactors = 2, option = "first")
Alpha:                 0.45 
G.6:                   0.55 
Omega Hierarchical:    0.56 
Omega H asymptotic:    0.89 
Omega Total            0.63 

Schmid Leiman Factor loadings greater than  0.2 
       g  F1*  F2*   h2   u2
V1            0.61 0.37 0.63
V2            0.63 0.40 0.60
V3  0.61           0.37 0.63
V4  0.64           0.41 0.59
V5-           0.58 0.34 0.66
V6            0.64 0.42 0.58
V7- 0.60           0.36 0.64
V8- 0.66           0.43 0.57

With eigenvalues of:
  g F1* F2* 
1.6 0.0 1.5 

general/max  1.04   max/min =   6.032862e+15
The degrees of freedom for the model is 13  and the fit was  0.06 
The number of observations was  500  with Chi Square =  31.63  with prob <  0.0027 

Measures of factor score adequacy              g F1* F2*
 Correlation of scores with factors            0.85 0 1.49
Multiple R square of scores with factors       0.73 0 2.21
Minimum correlation of factor score estimates  0.46 -1 3.41 
> #apply omega to analyze 6 mental ability tests 
> data(ability.cov)   #has a covariance matrix
> if(require(Rgraphviz)) {omega(ability.cov$cov)} else {omega(ability.cov$cov,plot=FALSE)} 
Omega 
Call: omega(m = ability.cov$cov)
Alpha:                 0.8 
G.6:                   0.83 
Omega Hierarchical:    0.57 
Omega H asymptotic:    0.64 
Omega Total            0.89 

Schmid Leiman Factor loadings greater than  0.2 
           g  F1*  F2*  F3*   h2   u2
general 0.60 0.35      0.20 0.56 0.44
picture 0.57      0.67      0.78 0.22
blocks  0.68      0.27 0.36 0.67 0.33
maze    0.45           0.46 0.42 0.58
reading 0.47 0.86           0.96     
vocab   0.45 0.68           0.66 0.34

With eigenvalues of:
   g  F1*  F2*  F3* 
1.78 1.32 0.57 0.38 

general/max  1.35   max/min =   3.44
The degrees of freedom for the model is 0  and the fit was  0 

Measures of factor score adequacy              g F1* F2* F3*
 Correlation of scores with factors            0.78 0.91 0.7 0.55
Multiple R square of scores with factors       0.61 0.82 0.5 0.3
Minimum correlation of factor score estimates  0.22 0.65 -0.01 -0.4 
> 
> 
> 
> cleanEx(); nameEx("omega.graph")
> ### * omega.graph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: omega.graph
> ### Title: Graph hierarchical factor structures
> ### Aliases: omega.graph omega.sem
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> #24 mental tests from Holzinger-Swineford-Harman
> if(require(GPArotation) ) {om24 <- omega(Harman74.cor$cov,4) } #run omega
Loading required package: GPArotation
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
> if(require(Rgraphviz) ){om24pn <- omega.graph(om24,sl=FALSE)} #show the structure
> #
> #example hierarchical structure from Jensen and Weng
> if(require(GPArotation) ) {jen.omega <- omega(make.hierarchical())}
Loading required package: MASS
> if(require(Rgraphviz) ) {om.jen <- omega.graph(jen.omega,sl=FALSE) }
> 
> 
> 
> 
> cleanEx(); nameEx("p.rep")
> ### * p.rep
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: p.rep
> ### Title: Find the probability of replication for an F, t, or r and
> ###   estimate effect size
> ### Aliases: p.rep p.rep.f p.rep.t p.rep.r
> ### Keywords: models univar
> 
> ### ** Examples
> 
> 
> p.rep(.05)  #probability of replicating a result if the original study had a  p = .05
[1] 0.877603
> p.rep.f(9.0,98)  #probability of replicating  a result with F = 9.0 with 98 df
$p.rep
[1] 0.9720729

$dprime
[1] 0.6060915

$prob
[1] 0.003423264

$r.equiv
[1] 0.2900209

> p.rep.r(.4,50)    #probability of replicating a result if r =.4 with n = 50
$p.rep
[1] 0.9818287

$dprime
[1] 0.8728716

$prob
[1] 0.003075876

> p.rep.t(3,98)   #probability of replicating a result if t = 3 with df =98
$p.rep
[1] 0.9720729

$dprime
[1] 0.6060915

$prob
[1] 0.003423264

$r.equiv
[1] 0.2900209

> p.rep.t(2.14,84,14) #effect of equal sample sizes (see Rosnow et al.)
$p.rep
[1] 0.9002023

$dprime
[1] 0.6054045

$prob
[1] 0.03483740

$r.equiv
[1] 0.2112921

> 
> 
> 
> 
> cleanEx(); nameEx("paired.r")
> ### * paired.r
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: paired.r
> ### Title: Test the difference between (un)paired correlations
> ### Aliases: paired.r
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> paired.r(.5,.3, .4, 100) #dependent correlations
$test
[1] "test of difference between two correlated  correlations"

$t
[1] 2.064994

$p
[1] 0.0415641

> paired.r(.5,.3,NULL,100) #independent correlations same sample size
$test
[1] "test of difference between two independent correlations"

$z
[1] 1.66992

$p
[1] 0.09493519

> paired.r(.5,.3,NULL, 100, 64) # independent correlations, different sample sizes
$test
[1] "test of difference between two independent correlations"

$z
[1] 1.467395

$p
[1] 0.1422686

> 
> 
> 
> cleanEx(); nameEx("pairs.panels")
> ### * pairs.panels
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pairs.panels
> ### Title: SPLOM, histograms and correlations for a data matrix
> ### Aliases: pairs.panels panel.cor panel.cor.scale panel.hist panel.lm
> ###   panel.lm.ellipse panel.hist.density panel.ellipse p.ellipse
> ###   panel.smoothie
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> 
> pairs.panels(attitude)   #see the graphics window
> data(peas)
> pairs.panels(peas,lm=TRUE,xlim=c(14,22),ylim=c(14,22))
> 
> 
> 
> 
> cleanEx(); nameEx("partial.r")
> ### * partial.r
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: partial.r
> ### Title: Find the partial correlations for a set (x) of variables with
> ###   set (y) removed.
> ### Aliases: partial.r
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> jen <- make.hierarchical()    #make up a correlation matrix 
Loading required package: MASS
> round(jen[1:5,1:5],2)
     V1   V2   V3   V4   V5
V1 1.00 0.56 0.48 0.40 0.35
V2 0.56 1.00 0.42 0.35 0.30
V3 0.48 0.42 1.00 0.30 0.26
V4 0.40 0.35 0.30 1.00 0.42
V5 0.35 0.30 0.26 0.42 1.00
> par.r <- partial.r(jen,c(1,3,5),c(2,4))
> par.r
     V1   V3   V5
V1 1.00 0.29 0.14
V3 0.29 1.00 0.10
V5 0.14 0.10 1.00
> 
> 
> 
> cleanEx(); nameEx("peas")
> ### * peas
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: peas
> ### Title: Galton`s Peas
> ### Aliases: peas
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(peas)
> pairs.panels(peas,lm=TRUE,xlim=c(14,22),ylim=c(14,22))
> describe(peas)
       var   n  mean   sd median trimmed  mad   min   max range skew kurtosis
parent   1 700 18.00 2.00  18.00   18.00 2.97 15.00 21.00   6.0 0.00    -1.25
child    2 700 16.29 1.98  16.07   16.14 2.55 13.77 22.67   8.9 0.49    -0.64
         se
parent 0.08
child  0.07
> pairs.panels(peas)
> 
> 
> 
> cleanEx(); nameEx("phi")
> ### * phi
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: phi
> ### Title: Find the phi coefficient of correlation between two dichotomous
> ###   variables
> ### Aliases: phi
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> phi(c(30,20,20,30))
[1] 0.2
> phi(c(40,10,10,40))
[1] 0.6
> x <- matrix(c(40,5,20,20),ncol=2)
> phi(x)
[1] 0.43
> 
> 
> 
> 
> cleanEx(); nameEx("phi.demo")
> ### * phi.demo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: phi.demo
> ### Title: A simple demonstration of the Pearson, phi, and polychoric
> ###   corelation
> ### Aliases: phi.demo
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> if(require(polycor)) {demo <- phi.demo() #compare the phi (lower off diagonal and polychoric correlations (upper off diagonal)
+ #show the result from poly.mat
+ round(demo$tetrachoric,2)
+ #show the result from phi2poly
+ #tetrachorics above the diagonal, phi below the diagonal 
+ round(demo$phis,2) }
Loading required package: polycor
Loading required package: mvtnorm
Loading required package: sfsmisc
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in log(P) : NaNs produced
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in hetcor.data.frame(dframe, ML = ML, std.err = std.err, bins = bins,  :
  the correlation matrix has been adjusted to make it positive-definite
         latent observed   X1   X2   X3   X4   X5
latent     1.00     0.60 0.26 0.39 0.47 0.39 0.26
observed   0.60     1.00 0.39 0.67 0.79 0.65 0.41
X1         0.26     0.39 1.00 0.99 0.98 0.94 0.89
X2         0.39     0.67 0.38 1.00 0.99 0.97 0.94
X3         0.47     0.79 0.17 0.44 1.00 0.99 0.97
X4         0.39     0.65 0.07 0.19 0.42 1.00 0.99
X5         0.26     0.41 0.03 0.07 0.17 0.40 1.00
> 
> 
> 
> cleanEx(); nameEx("phi2poly")
> ### * phi2poly
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: phi2poly
> ### Title: Convert a phi coefficient to a polychoric correlation
> ### Aliases: phi2poly
> ### Keywords: models models
> 
> ### ** Examples
> 
> #phi2poly(.3,.5,.5)
> #phi2poly(.3,.3,.7)
> 
> 
> 
> cleanEx(); nameEx("plot.psych")
> ### * plot.psych
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.psych
> ### Title: Plotting functions for the psych package of class ``psych"
> ### Aliases: plot.psych
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> test.data <- Harman74.cor$cov
> f4 <- factor.pa(test.data,4)
> plot(f4)
Use ICLUST.graph to see the  hierarchical structure
> 
> 
> 
> cleanEx(); nameEx("polar")
> ### * polar
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: polar
> ### Title: Convert Cartesian factor loadings into polar coordinates
> ### Aliases: polar
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> 
> circ.data <- circ.sim(24,500)
> circ.fa <- factor.pa(circ.data,2)
> circ.polar <- round(polar(circ.fa),2)
> circ.polar
    theta21 vl21
v8     2.74 0.40
v7    13.65 0.38
v6    30.64 0.41
v5    53.44 0.33
v4    66.25 0.38
v3    81.41 0.45
v2    96.55 0.35
v1   113.81 0.38
v24  120.09 0.37
v23  135.48 0.39
v22  155.76 0.36
v21  166.39 0.38
v20  182.52 0.34
v19  197.30 0.38
v18  213.84 0.38
v17  228.19 0.39
v16  248.04 0.36
v15  265.02 0.34
v14  273.06 0.39
v13  296.81 0.36
v12  303.55 0.44
v11  322.82 0.47
v10  337.78 0.34
v9   346.76 0.37
> #compare to the graphic
> cluster.plot(circ.fa)
> 
> 
> 
> cleanEx(); nameEx("poly.mat")
> ### * poly.mat
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: poly.mat
> ### Title: Find polychoric correlations of item data
> ### Aliases: poly.mat
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> if(require(polycor)) {demo <- phi.demo()
+ round(demo$tetrachoric,2) #these should actually be all 1s but won't be
+ round(demo$phis,2)    #compare the phi (lower off diagonal and polychoric correlations (upper off diagonal)
+ } else {print("I am sorry, this demo requires the polycor package")}
Loading required package: polycor
Loading required package: mvtnorm
Loading required package: sfsmisc
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in log(P) : NaNs produced
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in hetcor.data.frame(dframe, ML = ML, std.err = std.err, bins = bins,  :
  the correlation matrix has been adjusted to make it positive-definite
         latent observed   X1   X2   X3   X4   X5
latent     1.00     0.60 0.26 0.39 0.47 0.39 0.26
observed   0.60     1.00 0.39 0.67 0.79 0.65 0.41
X1         0.26     0.39 1.00 0.99 0.98 0.94 0.89
X2         0.39     0.67 0.38 1.00 0.99 0.97 0.94
X3         0.47     0.79 0.17 0.44 1.00 0.99 0.97
X4         0.39     0.65 0.07 0.19 0.42 1.00 0.99
X5         0.26     0.41 0.03 0.07 0.17 0.40 1.00
> 
> 
> 
> cleanEx(); nameEx("polychor.matrix")
> ### * polychor.matrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: polychor.matrix
> ### Title: Phi or Yule coefficient matrix to polychoric coefficient matrix
> ### Aliases: polychor.matrix Yule2poly.matrix phi2poly.matrix
> ###   Yule2phi.matrix
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> demo <- phi.demo()
Loading required package: polycor
Loading required package: mvtnorm
Loading required package: sfsmisc
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in log(P) : NaNs produced
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in optimise(f, interval = c(-1, 1)) :
  NA/Inf replaced by maximum positive value
Warning in hetcor.data.frame(dframe, ML = ML, std.err = std.err, bins = bins,  :
  the correlation matrix has been adjusted to make it positive-definite
> round(demo$tetrachoric,2) #these should actually be all 1s but won't be
     X1   X2   X3   X4   X5
X1 1.00 0.98 0.97 0.95 0.97
X2 0.98 1.00 0.99 0.96 0.95
X3 0.97 0.99 1.00 0.99 0.97
X4 0.95 0.96 0.99 1.00 0.98
X5 0.97 0.95 0.97 0.98 1.00
> round(demo$phis,2)    #compare the phi (lower off diagonal and polychoric correlations (upper off diagonal)
         latent observed   X1   X2   X3   X4   X5
latent     1.00     0.60 0.26 0.39 0.47 0.39 0.26
observed   0.60     1.00 0.39 0.67 0.79 0.65 0.41
X1         0.26     0.39 1.00 0.99 0.98 0.94 0.89
X2         0.39     0.67 0.38 1.00 0.99 0.97 0.94
X3         0.47     0.79 0.17 0.44 1.00 0.99 0.97
X4         0.39     0.65 0.07 0.19 0.42 1.00 0.99
X5         0.26     0.41 0.03 0.07 0.17 0.40 1.00
> 
> 
> 
> 
> cleanEx(); nameEx("principal")
> ### * principal
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: principal
> ### Title: Principal components analysis
> ### Aliases: principal
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> #Four principal components of the Harmon 24 variable problem
> #compare to a four factor principal axes solution using factor.congruence
> pc <- principal(Harman74.cor$cov,4,rotate="varimax")
> pa <- factor.pa(Harman74.cor$cov,4,rotate="varimax")
> round(factor.congruence(pc,pa),2)
     PA1  PA3  PA2  PA4
PC1 1.00 0.61 0.46 0.54
PC3 0.54 0.99 0.44 0.54
PC2 0.44 0.50 1.00 0.55
PC4 0.47 0.53 0.49 0.99
> 
> 
> 
> 
> cleanEx(); nameEx("print.psych")
> ### * print.psych
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: print.psych
> ### Title: Print and summary functions for the psych class
> ### Aliases: print.psych summary.psych
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(bfi)
>  keys.list <- list(agree=c(-1,2:5),conscientious=c(6:8,-9,-10),extraversion=c(-11,-12,13:15),neuroticism=c(16:20),openness = c(21,-22,23,24,-25))
>  keys <- make.keys(25,keys.list,item.labels=colnames(bfi))
>  scores <- score.items(keys,bfi,short=TRUE)
>  scores
Call: score.items(keys = keys, items = bfi, short = TRUE)

(Unstandardized) Alpha:
      agree conscientious extraversion neuroticism openness
alpha  0.65          0.74         0.77         0.8     0.47

Average item correlation:
          agree conscientious extraversion neuroticism openness
average.r  0.27          0.36         0.41        0.46     0.15

 Guttman 6* reliability: 
         agree conscientious extraversion neuroticism openness
Lambda.6  0.66          0.74         0.78        0.82     0.53

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
              agree conscientious extraversion neuroticism openness
agree          0.65          0.36         0.62       -0.30     0.33
conscientious  0.25          0.74         0.28       -0.22     0.24
extraversion   0.44          0.21         0.77       -0.24     0.57
neuroticism   -0.22         -0.17        -0.19        0.81     0.11
openness       0.18          0.14         0.34        0.07     0.47

Item by scale correlations:
 corrected for item overlap and scale reliability
   agree conscientious extraversion neuroticism openness
A1 -0.40         -0.12        -0.07        0.16    -0.16
A2  0.63          0.20         0.41       -0.09     0.26
A3  0.63          0.20         0.41       -0.13     0.30
A4  0.39          0.21         0.22       -0.17    -0.05
A5  0.61          0.19         0.54       -0.23     0.17
C1  0.12          0.59         0.11       -0.03     0.23
C2  0.19          0.56         0.08        0.07     0.14
C3  0.24          0.61         0.23       -0.03     0.15
C4 -0.29         -0.70        -0.22        0.30    -0.17
C5 -0.22         -0.57        -0.20        0.32    -0.01
E1 -0.36         -0.07        -0.65        0.09    -0.42
E2 -0.35         -0.21        -0.69        0.28    -0.23
E3  0.44          0.19         0.61       -0.09     0.50
E4  0.53          0.12         0.68       -0.23     0.10
E5  0.27          0.34         0.58       -0.05     0.49
N1 -0.26         -0.17        -0.06        0.77     0.13
N2 -0.26         -0.11        -0.06        0.74     0.09
N3 -0.19         -0.16        -0.10        0.71     0.10
N4 -0.24         -0.27        -0.37        0.63     0.08
N5 -0.06         -0.04        -0.23        0.54    -0.03
O1  0.21          0.17         0.30       -0.01     0.54
O2  0.13          0.13        -0.22       -0.06    -0.09
O3  0.29          0.21         0.43        0.01     0.63
O4  0.17          0.10         0.00        0.22     0.42
O5 -0.12         -0.14        -0.15        0.06    -0.42
>  summary(scores)
Call: score.items(keys = keys, items = bfi, short = TRUE)

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, (unstandardized) alpha on the diagonal 
 corrected correlations above the diagonal:
              agree conscientious extraversion neuroticism openness
agree          0.65          0.36         0.62       -0.30     0.33
conscientious  0.25          0.74         0.28       -0.22     0.24
extraversion   0.44          0.21         0.77       -0.24     0.57
neuroticism   -0.22         -0.17        -0.19        0.81     0.11
openness       0.18          0.14         0.34        0.07     0.47
> 
> 
> 
> cleanEx(); nameEx("r.test")
> ### * r.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: r.test
> ### Title: Tests of significance for correlations
> ### Aliases: r.test
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> n <- 30 
> r <- seq(0,.9,.1) 
> rc <- matrix(r.con(r,n),ncol=2) 
> test <- r.test(n,r)
> r.rc <- data.frame(r=r,z=fisherz(r),lower=rc[,1],upper=rc[,2],t=test$t,p=test$p) 
> round(r.rc,2) 
     r    z lower upper     t    p
1  0.0 0.00 -0.36  0.36  0.00 1.00
2  0.1 0.10 -0.27  0.44  0.53 0.60
3  0.2 0.20 -0.17  0.52  1.08 0.29
4  0.3 0.31 -0.07  0.60  1.66 0.11
5  0.4 0.42  0.05  0.66  2.31 0.03
6  0.5 0.55  0.17  0.73  3.06 0.00
7  0.6 0.69  0.31  0.79  3.97 0.00
8  0.7 0.87  0.45  0.85  5.19 0.00
9  0.8 1.10  0.62  0.90  7.06 0.00
10 0.9 1.47  0.80  0.95 10.93 0.00
> 
> r.test(50,r)
Correlation tests 
Call:r.test(n = 50, r12 = r)
Test of significance of a  correlation 
 t value 0 0.7 1.41 2.18 3.02 4 5.2 6.79 9.24 14.3    with probability < 1 0.49 0.16 0.034 0.004 0.00022 4.1e-06 1.5e-08 3.2e-12 0
 and confidence interval  -0.28 -0.18 -0.08 0.02 0.14 0.26 0.39 0.52 0.67 0.83 0.28 0.37 0.45 0.53 0.61 0.68 0.75 0.82 0.88 0.94> r.test(30,.4,.6)       #test the difference between two independent correlations
Correlation tests 
Call:r.test(n = 30, r12 = 0.4, r34 = 0.6)
Test of difference between two independent correlations 
 z value 0.99    with probability  0.32> r.test(103,.4,.5,.1)   #Steiger case A
Correlation tests 
Call:r.test(n = 103, r12 = 0.4, r34 = 0.5, r23 = 0.1)
Test of difference between two correlated  correlations 
 t value -0.89    with probability < 0.37> r.test(103,.5,.6,.7,.5,.5,.8)  #steiger Case B 
Correlation tests 
Call:r.test(n = 103, r12 = 0.5, r34 = 0.6, r23 = 0.7, r13 = 0.5, r14 = 0.5, 
    r24 = 0.8)
Test of difference between two dependent correlations 
 z value -1.2    with probability  0.23> 
> 
> 
> 
> cleanEx(); nameEx("read.clipboard")
> ### * read.clipboard
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: read.clipboard
> ### Title: shortcut for reading from the clipboard
> ### Aliases: read.clipboard read.clipboard.csv read.clipboard.tab
> ###   read.clipboard.lower read.clipboard.upper
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> #my.data <- read.clipboad()
> #my.data <- read.clipboard.csv()
> #my.data <- read.clipboad(header=FALSE)
> 
> 
> 
> cleanEx(); nameEx("rescale")
> ### * rescale
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rescale
> ### Title: Function to convert scores to ``conventional " metrics
> ### Aliases: rescale
> ### Keywords: multivariate models univar
> 
> ### ** Examples
> 
> T <- rescale(attitude,50,10) #all put on same scale
> describe(T)
           var  n mean sd median trimmed   mad   min   max range  skew kurtosis
rating       1 30   50 10  50.71   50.47  8.53 29.76 66.73 36.97 -0.36    -0.77
complaints   2 30   50 10  48.80   50.36 11.14 27.77 67.57 39.81 -0.22    -0.68
privileges   3 30   50 10  48.67   49.69  8.48 31.09 74.41 43.32  0.38    -0.41
learning     4 30   50 10  50.11   50.18 12.63 30.94 65.88 34.93 -0.05    -1.22
raises       5 30   50 10  48.91   49.87 10.69 29.19 72.47 43.28  0.20    -0.60
critical     6 30   50 10  52.76   51.08  7.49 23.96 67.42 43.46 -0.87     0.17
advance      7 30   50 10  48.12   48.93  8.65 32.57 78.25 45.68  0.85     0.47
             se
rating     1.83
complaints 1.83
privileges 1.83
learning   1.83
raises     1.83
critical   1.83
advance    1.83
> T1 <- rescale(attitude,seq(0,300,50),seq(10,70,10)) #different means and sigmas
> describe(T1)
           var  n mean sd median trimmed   mad    min    max  range  skew
rating       1 30    0 10   0.71    0.47  8.53 -20.24  16.73  36.97 -0.36
complaints   2 30   50 20  47.60   50.73 22.27   5.54  85.15  79.61 -0.22
privileges   3 30  100 30  96.00   99.06 25.45  43.28 173.23 129.95  0.38
learning     4 30  150 40 150.45  150.74 50.53  73.77 213.50 139.73 -0.05
raises       5 30  200 50 194.55  199.36 53.47  95.97 312.37 216.40  0.20
critical     6 30  250 60 266.57  256.47 44.95  93.76 354.50 260.74 -0.87
advance      7 30  300 70 286.85  292.52 60.52 177.99 497.76 319.77  0.85
           kurtosis    se
rating        -0.77  1.83
complaints    -0.68  3.65
privileges    -0.41  5.48
learning      -1.22  7.30
raises        -0.60  9.13
critical       0.17 10.95
advance        0.47 12.78
> 
> 
> 
> cleanEx(); nameEx("sat.act")
> ### * sat.act
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sat.act
> ### Title: 3 Measures of ability: SATV, SATQ, ACT
> ### Aliases: sat.act
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(sat.act)
> describe(sat.act)
          var   n   mean     sd median trimmed    mad min max range  skew
gender      1 700   1.65   0.48      2    1.68   0.00   1   2     1 -0.61
education   2 700   3.16   1.43      3    3.31   1.48   0   5     5 -0.68
age         3 700  25.59   9.50     22   23.86   5.93  13  65    52  1.64
ACT         4 700  28.55   4.82     29   28.84   4.45   3  36    33 -0.66
SATV        5 700 612.23 112.90    620  619.45 118.61 200 800   600 -0.64
SATQ        6 687 610.22 115.64    620  617.25 118.61 200 800   600 -0.59
          kurtosis   se
gender       -1.62 0.02
education    -0.07 0.05
age           2.42 0.36
ACT           0.53 0.18
SATV          0.33 4.27
SATQ         -0.02 4.41
> pairs.panels(sat.act)
> 
> 
> 
> cleanEx(); nameEx("schmid")
> ### * schmid
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: schmid
> ### Title: Apply the Schmid Leiman transformation to a correlation matrix
> ### Aliases: schmid
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> jen <- sim.hierarchical()  #create a hierarchical demo
Loading required package: MASS
> o.jen <- schmid(jen,digits=2)  #use the oblimin rotation
Loading required package: GPArotation
> p.jen <- schmid(jen,rotate="promax") #use the promax rotation
> 
> 
> 
> cleanEx(); nameEx("score.alpha")
> ### * score.alpha
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: score.alpha
> ### Title: Score scales and find Cronbach's alpha as well as associated
> ###   statistics
> ### Aliases: score.alpha
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> y <- attitude     #from the datasets package
> keys <- matrix(c(rep(1,7),rep(1,4),rep(0,7),rep(-1,3)),ncol=3)
> labels <- c("first","second","third")
> x <- score.alpha(keys,y,labels)
> 
> 
> 
> 
> cleanEx(); nameEx("score.items")
> ### * score.items
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: score.items
> ### Title: Score item composite scales and find Cronbach's alpha, Guttman
> ###   lambda 6 and item whole correlations
> ### Aliases: score.items
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> 
> #see  the example including the bfi data set
> data(bfi)
>  keys.list <- list(agree=c(-1,2:5),conscientious=c(6:8,-9,-10),extraversion=c(-11,-12,13:15),neuroticism=c(16:20),openness = c(21,-22,23,24,-25))
>  keys <- make.keys(25,keys.list,item.labels=colnames(bfi))
>  scores <- score.items(keys,bfi)
>  scores
Call: score.items(keys = keys, items = bfi)

(Unstandardized) Alpha:
      agree conscientious extraversion neuroticism openness
alpha  0.65          0.74         0.77        0.81     0.47

Average item correlation:
          agree conscientious extraversion neuroticism openness
average.r  0.27          0.36         0.41        0.46     0.15

 Guttman 6* reliability: 
         agree conscientious extraversion neuroticism openness
Lambda.6  0.66          0.74         0.78        0.82     0.53

Scale intercorrelations corrected for attenuation 
 raw correlations below the diagonal, alpha on the diagonal 
 corrected correlations above the diagonal:
              agree conscientious extraversion neuroticism openness
agree          0.65          0.36         0.62       -0.30     0.33
conscientious  0.25          0.74         0.28       -0.22     0.24
extraversion   0.44          0.21         0.77       -0.24     0.57
neuroticism   -0.22         -0.17        -0.19        0.81     0.11
openness       0.18          0.14         0.34        0.07     0.47

Item by scale correlations:
 corrected for item overlap and scale reliability
   agree conscientious extraversion neuroticism openness
A1 -0.40         -0.12        -0.07        0.16    -0.16
A2  0.63          0.20         0.41       -0.09     0.26
A3  0.63          0.20         0.41       -0.13     0.30
A4  0.39          0.21         0.22       -0.17    -0.05
A5  0.61          0.19         0.54       -0.23     0.17
C1  0.12          0.59         0.11       -0.03     0.23
C2  0.19          0.56         0.08        0.07     0.14
C3  0.24          0.61         0.23       -0.03     0.15
C4 -0.29         -0.70        -0.22        0.30    -0.17
C5 -0.22         -0.57        -0.20        0.32    -0.01
E1 -0.36         -0.07        -0.65        0.09    -0.42
E2 -0.35         -0.21        -0.69        0.28    -0.23
E3  0.44          0.19         0.61       -0.09     0.50
E4  0.53          0.12         0.68       -0.23     0.10
E5  0.27          0.34         0.58       -0.05     0.49
N1 -0.26         -0.17        -0.06        0.77     0.13
N2 -0.26         -0.11        -0.06        0.74     0.09
N3 -0.19         -0.16        -0.10        0.71     0.10
N4 -0.24         -0.27        -0.37        0.63     0.08
N5 -0.06         -0.04        -0.23        0.54    -0.03
O1  0.21          0.17         0.30       -0.01     0.54
O2  0.13          0.13        -0.22       -0.06    -0.09
O3  0.29          0.21         0.43        0.01     0.63
O4  0.17          0.10         0.00        0.22     0.42
O5 -0.12         -0.14        -0.15        0.06    -0.42
> 
> 
> 
> 
> cleanEx(); nameEx("score.multiple.choice")
> ### * score.multiple.choice
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: score.multiple.choice
> ### Title: Score multiple choice items and provide basic test statistics
> ### Aliases: score.multiple.choice
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> data(iqitems)
> iq.keys <- c(4,4,3,1,4,3,2,3,1,4,1,3,4,3)
> score.multiple.choice(iq.keys,iqitems)
$item.stats
     key    0    1    2    3    4    5    6    r    n mean   sd  skew kurtosis
iq1    4 0.04 0.01 0.03 0.09 0.80 0.02 0.01 0.59 1000 0.80 0.40 -1.51     0.27
iq8    4 0.03 0.10 0.01 0.02 0.80 0.01 0.04 0.39 1000 0.80 0.40 -1.49     0.22
iq10   3 0.10 0.22 0.09 0.37 0.04 0.13 0.04 0.35 1000 0.37 0.48  0.53    -1.72
iq15   1 0.03 0.65 0.16 0.15 0.00 0.00 0.00 0.35 1000 0.65 0.48 -0.63    -1.60
iq20   4 0.03 0.02 0.03 0.03 0.85 0.02 0.01 0.42 1000 0.85 0.35 -2.00     2.01
iq44   3 0.03 0.10 0.06 0.64 0.02 0.14 0.01 0.42 1000 0.64 0.48 -0.61    -1.64
iq47   2 0.04 0.08 0.59 0.06 0.11 0.07 0.05 0.51 1000 0.59 0.49 -0.35    -1.88
iq2    3 0.07 0.08 0.31 0.32 0.15 0.05 0.02 0.26 1000 0.32 0.46  0.80    -1.37
iq11   1 0.04 0.87 0.03 0.01 0.01 0.01 0.04 0.54 1000 0.87 0.34 -2.15     2.61
iq16   4 0.05 0.05 0.08 0.07 0.74 0.01 0.00 0.56 1000 0.74 0.44 -1.11    -0.77
iq32   1 0.04 0.54 0.02 0.14 0.10 0.04 0.12 0.50 1000 0.54 0.50 -0.17    -1.97
iq37   3 0.07 0.10 0.09 0.26 0.13 0.02 0.34 0.23 1000 0.26 0.44  1.12    -0.74
iq43   4 0.04 0.07 0.04 0.02 0.78 0.03 0.00 0.50 1000 0.78 0.41 -1.35    -0.18
iq49   3 0.06 0.27 0.09 0.32 0.14 0.08 0.05 0.28 1000 0.32 0.47  0.79    -1.38
       se
iq1  0.01
iq8  0.01
iq10 0.02
iq15 0.02
iq20 0.01
iq44 0.02
iq47 0.02
iq2  0.01
iq11 0.01
iq16 0.01
iq32 0.02
iq37 0.01
iq43 0.01
iq49 0.01

$alpha
         Averages
Averages     0.63

$av.r
         Averages
Averages     0.11

> #just convert the items to true or false 
> iq.tf <- score.multiple.choice(iq.keys,iqitems,score=FALSE)
> describe(iq.tf)  #compare to previous results
     var    n mean   sd median trimmed mad min max range  skew kurtosis   se
iq1    1 1000 0.80 0.40      1    0.88   0   0   1     1 -1.51     0.27 0.01
iq8    2 1000 0.80 0.40      1    0.87   0   0   1     1 -1.49     0.22 0.01
iq10   3 1000 0.37 0.48      0    0.34   0   0   1     1  0.53    -1.72 0.02
iq15   4 1000 0.65 0.48      1    0.69   0   0   1     1 -0.63    -1.60 0.02
iq20   5 1000 0.85 0.35      1    0.94   0   0   1     1 -2.00     2.01 0.01
iq44   6 1000 0.64 0.48      1    0.68   0   0   1     1 -0.61    -1.64 0.02
iq47   7 1000 0.59 0.49      1    0.61   0   0   1     1 -0.35    -1.88 0.02
iq2    8 1000 0.32 0.46      0    0.27   0   0   1     1  0.80    -1.37 0.01
iq11   9 1000 0.87 0.34      1    0.96   0   0   1     1 -2.15     2.61 0.01
iq16  10 1000 0.74 0.44      1    0.80   0   0   1     1 -1.11    -0.77 0.01
iq32  11 1000 0.54 0.50      1    0.55   0   0   1     1 -0.17    -1.97 0.02
iq37  12 1000 0.26 0.44      0    0.19   0   0   1     1  1.12    -0.74 0.01
iq43  13 1000 0.78 0.41      1    0.85   0   0   1     1 -1.35    -0.18 0.01
iq49  14 1000 0.32 0.47      0    0.27   0   0   1     1  0.79    -1.38 0.01
> 
> 
> 
> 
> cleanEx(); nameEx("sim")
> ### * sim
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim
> ### Title: Functions to simulate psychological/psychometric data.
> ### Aliases: sim sim.simplex
> ### Keywords: multivariate datagen
> 
> ### ** Examples
> 
> simplex <- sim()
Loading required package: MASS
> round(simplex$model,2)
      V1   V2   V3   V4   V5   V6   V7   V8   V9  V10  V11  V12
V1  1.00 0.56 0.48 0.51 0.45 0.38 0.41 0.36 0.31 0.33 0.29 0.25
V2  0.56 1.00 0.42 0.45 0.39 0.34 0.36 0.31 0.27 0.29 0.25 0.22
V3  0.48 0.42 1.00 0.38 0.34 0.29 0.31 0.27 0.23 0.25 0.22 0.18
V4  0.51 0.45 0.38 1.00 0.56 0.48 0.51 0.45 0.38 0.41 0.36 0.31
V5  0.45 0.39 0.34 0.56 1.00 0.42 0.45 0.39 0.34 0.36 0.31 0.27
V6  0.38 0.34 0.29 0.48 0.42 1.00 0.38 0.34 0.29 0.31 0.27 0.23
V7  0.41 0.36 0.31 0.51 0.45 0.38 1.00 0.56 0.48 0.51 0.45 0.38
V8  0.36 0.31 0.27 0.45 0.39 0.34 0.56 1.00 0.42 0.45 0.39 0.34
V9  0.31 0.27 0.23 0.38 0.34 0.29 0.48 0.42 1.00 0.38 0.34 0.29
V10 0.33 0.29 0.25 0.41 0.36 0.31 0.51 0.45 0.38 1.00 0.56 0.48
V11 0.29 0.25 0.22 0.36 0.31 0.27 0.45 0.39 0.34 0.56 1.00 0.42
V12 0.25 0.22 0.18 0.31 0.27 0.23 0.38 0.34 0.29 0.48 0.42 1.00
> 
> congeneric <- sim.congeneric()
> round(congeneric,2)
     V1   V2   V3   V4
V1 1.00 0.56 0.48 0.40
V2 0.56 1.00 0.42 0.35
V3 0.48 0.42 1.00 0.30
V4 0.40 0.35 0.30 1.00
> R <- sim.hierarchical()
> R
       V1     V2     V3     V4     V5    V6     V7     V8     V9
V1 1.0000 0.5600 0.4800 0.4032 0.3456 0.288 0.3024 0.2520 0.2016
V2 0.5600 1.0000 0.4200 0.3528 0.3024 0.252 0.2646 0.2205 0.1764
V3 0.4800 0.4200 1.0000 0.3024 0.2592 0.216 0.2268 0.1890 0.1512
V4 0.4032 0.3528 0.3024 1.0000 0.4200 0.350 0.2352 0.1960 0.1568
V5 0.3456 0.3024 0.2592 0.4200 1.0000 0.300 0.2016 0.1680 0.1344
V6 0.2880 0.2520 0.2160 0.3500 0.3000 1.000 0.1680 0.1400 0.1120
V7 0.3024 0.2646 0.2268 0.2352 0.2016 0.168 1.0000 0.3000 0.2400
V8 0.2520 0.2205 0.1890 0.1960 0.1680 0.140 0.3000 1.0000 0.2000
V9 0.2016 0.1764 0.1512 0.1568 0.1344 0.112 0.2400 0.2000 1.0000
> fx <- matrix(c(.9,.8,.7,rep(0,6),c(.8,.7,.6)),ncol=2)
> fy <- c(.6,.5,.4)
> Phi <- matrix(c(1,0,.5,0,1,.4,0,0,0),ncol=3)
> print(sim.structure(fx,Phi,fy,),digits=2)  
Call: NULL

 $model (Population correlation matrix) 
     V1   V2   V3   V4   V5    V6 V7 V8 V9
V1 1.00 0.72 0.63 0.00 0.00 0.000  0  0  0
V2 0.72 1.00 0.56 0.00 0.00 0.000  0  0  0
V3 0.63 0.56 1.00 0.00 0.00 0.000  0  0  0
V4 0.00 0.00 0.00 1.00 0.56 0.480  0  0  0
V5 0.00 0.00 0.00 0.56 1.00 0.420  0  0  0
V6 0.00 0.00 0.00 0.48 0.42 1.000  0  0  0
V7 0.27 0.24 0.21 0.19 0.17 0.144  1  0  0
V8 0.23 0.20 0.17 0.16 0.14 0.120  0  1  0
V9 0.18 0.16 0.14 0.13 0.11 0.096  0  0  1

$reliability (population reliability) 
 Vx1  Vx2  Vx3  Vx4  Vx5  Vx6  Vy1  Vy2  Vy3 
0.81 0.64 0.49 0.64 0.49 0.36 0.36 0.25 0.16 
> cor.plot(R) #show it graphically
> 
> simp <- sim.simplex()
> #show the simplex structure using cor.plot
> cor.plot(simp)
> 
> 
> 
> cleanEx(); nameEx("sim.VSS")
> ### * sim.VSS
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.VSS
> ### Title: create VSS like data
> ### Aliases: sim.VSS VSS.simulate VSS.sim
> ### Keywords: multivariate models datagen
> 
> ### ** Examples
> 
> ## Not run: 
> ##D simulated <- sim.VSS(1000,20,4,.6)
> ##D vss <- VSS(simulated,rotate="varimax")
> ##D VSS.plot(vss)
> ## End(Not run)
> 
> 
> 
> 
> cleanEx(); nameEx("sim.anova")
> ### * sim.anova
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.anova
> ### Title: Simulate a 3 way balanced ANOVA or linear model, with or without
> ###   repeated measures.
> ### Aliases: sim.anova
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> set.seed(42)
> data.df <- sim.anova(es1=1,es2=.5,es13=1)  # one main effect and one interaction
> describe(data.df)
     var  n mean   sd median trimmed  mad   min  max range  skew kurtosis   se
IV1*   1 16 1.50 0.52    1.5    1.50 0.74  1.00 2.00   1.0  0.00    -2.12 0.13
IV2*   2 16 1.50 0.52    1.5    1.50 0.74  1.00 2.00   1.0  0.00    -2.12 0.13
IV3*   3 16 1.50 0.52    1.5    1.50 0.74  1.00 2.00   1.0  0.00    -2.12 0.13
DV     4 16 0.49 1.85    1.0    0.62 1.63 -3.78 3.03   6.8 -0.67    -0.45 0.46
> pairs.panels(data.df)   #show how the design variables are orthogonal
> #
> summary(lm(DV~IV1*IV2*IV3,data=data.df))

Call:
lm(formula = DV ~ IV1 * IV2 * IV3, data = data.df)

Residuals:
       Min         1Q     Median         3Q        Max 
-8.966e-01 -3.917e-01  2.498e-16  3.917e-01  8.966e-01 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1.1798     0.5803   2.033 0.076503 .  
IV11            -1.9469     0.8207  -2.372 0.045093 *  
IV21             0.1076     0.8207   0.131 0.898974    
IV31            -4.0620     0.8207  -4.949 0.001122 ** 
IV11:IV21        2.6342     1.1607   2.269 0.052933 .  
IV11:IV31        6.0582     1.1607   5.220 0.000803 ***
IV21:IV31        2.0421     1.1607   1.759 0.116556    
IV11:IV21:IV31  -3.3524     1.6415  -2.042 0.075397 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.8207 on 8 degrees of freedom
Multiple R-squared: 0.8952,	Adjusted R-squared: 0.8035 
F-statistic: 9.764 on 7 and 8 DF,  p-value: 0.002274 

> summary(aov(DV~IV1*IV2*IV3,data=data.df))
            Df  Sum Sq Mean Sq F value    Pr(>F)    
IV1          1  9.7491  9.7491 14.4733 0.0052041 ** 
IV2          1 10.3370 10.3370 15.3461 0.0044345 ** 
IV3          1  2.8896  2.8896  4.2898 0.0720970 .  
IV1:IV2      1  0.9176  0.9176  1.3623 0.2767507    
IV1:IV3      1 19.2020 19.2020 28.5069 0.0006948 ***
IV2:IV3      1  0.1339  0.1339  0.1987 0.6675723    
IV1:IV2:IV3  1  2.8097  2.8097  4.1712 0.0753968 .  
Residuals    8  5.3887  0.6736                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> set.seed(42)
> data.df <- sim.anova(es1=1,es2=.5,es13=1,center=FALSE)  # demonstrate the effect of not centering the data on the regression
> describe(data.df)
     var  n mean   sd median trimmed  mad   min  max range  skew kurtosis   se
IV1*   1 16 1.50 0.52    1.5    1.50 0.74  1.00 2.00   1.0  0.00    -2.12 0.13
IV2*   2 16 1.50 0.52    1.5    1.50 0.74  1.00 2.00   1.0  0.00    -2.12 0.13
IV3*   3 16 1.50 0.52    1.5    1.50 0.74  1.00 2.00   1.0  0.00    -2.12 0.13
DV     4 16 0.49 1.85    1.0    0.62 1.63 -3.78 3.03   6.8 -0.67    -0.45 0.46
> #
> summary(lm(DV~IV1*IV2*IV3,data=data.df)) #this one is incorrect, because the IVs are not centered

Call:
lm(formula = DV ~ IV1 * IV2 * IV3, data = data.df)

Residuals:
       Min         1Q     Median         3Q        Max 
-8.966e-01 -3.917e-01  2.498e-16  3.917e-01  8.966e-01 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1.1798     0.5803   2.033 0.076503 .  
IV11            -1.9469     0.8207  -2.372 0.045093 *  
IV21             0.1076     0.8207   0.131 0.898974    
IV31            -4.0620     0.8207  -4.949 0.001122 ** 
IV11:IV21        2.6342     1.1607   2.269 0.052933 .  
IV11:IV31        6.0582     1.1607   5.220 0.000803 ***
IV21:IV31        2.0421     1.1607   1.759 0.116556    
IV11:IV21:IV31  -3.3524     1.6415  -2.042 0.075397 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.8207 on 8 degrees of freedom
Multiple R-squared: 0.8952,	Adjusted R-squared: 0.8035 
F-statistic: 9.764 on 7 and 8 DF,  p-value: 0.002274 

> summary(aov(DV~IV1*IV2*IV3,data=data.df)) #compare with the lm model
            Df  Sum Sq Mean Sq F value    Pr(>F)    
IV1          1  9.7491  9.7491 14.4733 0.0052041 ** 
IV2          1 10.3370 10.3370 15.3461 0.0044345 ** 
IV3          1  2.8896  2.8896  4.2898 0.0720970 .  
IV1:IV2      1  0.9176  0.9176  1.3623 0.2767507    
IV1:IV3      1 19.2020 19.2020 28.5069 0.0006948 ***
IV2:IV3      1  0.1339  0.1339  0.1987 0.6675723    
IV1:IV2:IV3  1  2.8097  2.8097  4.1712 0.0753968 .  
Residuals    8  5.3887  0.6736                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> #now examine multiple levels and quadratic terms
> set.seed(42)
> data.df <- sim.anova(es1=1,es13=1,n2=3,n3=4,es22=1)
> summary(lm(DV~IV1*IV2*IV3,data=data.df))

Call:
lm(formula = DV ~ IV1 * IV2 * IV3, data = data.df)

Residuals:
       Min         1Q     Median         3Q        Max 
-2.002e+00 -3.397e-01  5.638e-16  3.397e-01  2.002e+00 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       3.4260     0.7844   4.367 0.000207 ***
IV11             -2.7790     1.1093  -2.505 0.019434 *  
IV20             -3.0489     1.1093  -2.748 0.011190 *  
IV21             -1.2009     1.1093  -1.083 0.289778    
IV3-1            -1.5254     1.1093  -1.375 0.181825    
IV31             -4.4713     1.1093  -4.031 0.000488 ***
IV33             -5.1016     1.1093  -4.599 0.000115 ***
IV11:IV20         1.5126     1.5689   0.964 0.344596    
IV11:IV21         1.3254     1.5689   0.845 0.406549    
IV11:IV3-1        3.2038     1.5689   2.042 0.052272 .  
IV11:IV31         6.1556     1.5689   3.924 0.000639 ***
IV11:IV33         8.5233     1.5689   5.433 1.40e-05 ***
IV20:IV3-1        2.1234     1.5689   1.353 0.188512    
IV21:IV3-1        1.1223     1.5689   0.715 0.481280    
IV20:IV31         1.3930     1.5689   0.888 0.383389    
IV21:IV31         2.2484     1.5689   1.433 0.164711    
IV20:IV33         1.5838     1.5689   1.010 0.322782    
IV21:IV33         1.5504     1.5689   0.988 0.332899    
IV11:IV20:IV3-1  -2.6968     2.2187  -1.215 0.236008    
IV11:IV21:IV3-1  -1.2671     2.2187  -0.571 0.573237    
IV11:IV20:IV31   -0.4246     2.2187  -0.191 0.849846    
IV11:IV21:IV31   -3.3169     2.2187  -1.495 0.147959    
IV11:IV20:IV33   -2.4872     2.2187  -1.121 0.273368    
IV11:IV21:IV33   -0.6422     2.2187  -0.289 0.774713    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.109 on 24 degrees of freedom
Multiple R-squared: 0.8651,	Adjusted R-squared: 0.7359 
F-statistic: 6.693 on 23 and 24 DF,  p-value: 8.495e-06 

> summary(aov(DV~IV1*IV2*IV3,data=data.df))
            Df Sum Sq Mean Sq F value    Pr(>F)    
IV1          1 36.115  36.115 29.3460 1.451e-05 ***
IV2          2 32.233  16.116 13.0958 0.0001429 ***
IV3          3 10.759   3.586  2.9142 0.0549549 .  
IV1:IV2      2  0.028   0.014  0.0113 0.9887276    
IV1:IV3      3 98.019  32.673 26.5493 8.488e-08 ***
IV2:IV3      6  3.707   0.618  0.5021 0.8005428    
IV1:IV2:IV3  6  8.581   1.430  1.1621 0.3586822    
Residuals   24 29.536   1.231                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> pairs.panels(data.df)
> #
> data.df <- sim.anova(es1=1,es2=-.5,within=c(-1,0,1),n=10)
> pairs.panels(data.df)
> 
> 
> 
> 
> cleanEx(); nameEx("sim.congeneric")
> ### * sim.congeneric
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.congeneric
> ### Title: Simulate a congeneric data set
> ### Aliases: congeneric.sim sim.congeneric make.congeneric
> ### Keywords: multivariate datagen
> 
> ### ** Examples
> 
> test <- sim.congeneric(c(.9,.8,.7,.6))   #just the population matrix
> test <- sim.congeneric(c(.9,.8,.7,.6),N=100)   # a sample correlation matrix
> test <- sim.congeneric(short=FALSE, N=100)
> round(cor(test$observed),2) # show  a congeneric correlation matrix
     V1   V2   V3   V4
V1 1.00 0.50 0.49 0.28
V2 0.50 1.00 0.37 0.29
V3 0.49 0.37 1.00 0.42
V4 0.28 0.29 0.42 1.00
> f1=factor.pa(test$observed,1,scores=TRUE)
> round(cor(f1$scores,test$latent),2)  #factor score estimates are correlated with but not equal to the factor scores
    theta   e1   e2   e3   e4
PA1  0.84 0.29 0.14 0.47 0.18
> 
> 
> 
> 
> 
> cleanEx(); nameEx("sim.hierarchical")
> ### * sim.hierarchical
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.hierarchical
> ### Title: Create a population or sample correlation matrix, perhaps with
> ###   hierarchical structure.
> ### Aliases: sim.hierarchical make.hierarchical
> ### Keywords: multivariate models datagen
> 
> ### ** Examples
> 
> 
> gload <-  gload<-matrix(c(.9,.8,.7),nrow=3)    # a higher order factor matrix
> fload <-matrix(c(                    #a lower order (oblique) factor matrix
+            .8,0,0,
+            .7,0,.0,
+            .6,0,.0,
+             0,.7,.0,
+             0,.6,.0,
+             0,.5,0,
+             0,0,.6,
+             0,0,.5,
+             0,0,.4),   ncol=3,byrow=TRUE)
>             
> jensen <- sim.hierarchical(gload,fload)    #the test set used by omega
Loading required package: MASS
> round(jensen,2)     
     V1   V2   V3   V4   V5   V6   V7   V8   V9
V1 1.00 0.56 0.48 0.40 0.35 0.29 0.30 0.25 0.20
V2 0.56 1.00 0.42 0.35 0.30 0.25 0.26 0.22 0.18
V3 0.48 0.42 1.00 0.30 0.26 0.22 0.23 0.19 0.15
V4 0.40 0.35 0.30 1.00 0.42 0.35 0.24 0.20 0.16
V5 0.35 0.30 0.26 0.42 1.00 0.30 0.20 0.17 0.13
V6 0.29 0.25 0.22 0.35 0.30 1.00 0.17 0.14 0.11
V7 0.30 0.26 0.23 0.24 0.20 0.17 1.00 0.30 0.24
V8 0.25 0.22 0.19 0.20 0.17 0.14 0.30 1.00 0.20
V9 0.20 0.18 0.15 0.16 0.13 0.11 0.24 0.20 1.00
> 
> fload <- matrix(c(c(c(.9,.8,.7,.6),rep(0,20)),c(c(.9,.8,.7,.6),rep(0,20)),c(c(.9,.8,.7,.6),rep(0,20)),c(c(c(.9,.8,.7,.6),rep(0,20)),c(.9,.8,.7,.6))),ncol=5)
> gload <- matrix(rep(0,5))
> five.factor <- sim.hierarchical(gload,fload,500,TRUE) #create sample data set
> 
> 
> 
> 
> cleanEx(); nameEx("sim.item")
> ### * sim.item
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.item
> ### Title: Generate simulated data structures for circumplex or simple
> ###   structure
> ### Aliases: item.sim sim.item sim.dichot item.dichot sim.circ circ.sim
> ### Keywords: multivariate datagen
> 
> ### ** Examples
> 
> 
> round(cor(circ.sim(nvar=8,nsub=200)),2)
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,]  1.00  0.26 -0.04 -0.23 -0.21 -0.22  0.11  0.28
[2,]  0.26  1.00  0.21  0.00 -0.14 -0.35 -0.20  0.08
[3,] -0.04  0.21  1.00  0.35  0.08 -0.26 -0.42 -0.27
[4,] -0.23  0.00  0.35  1.00  0.23 -0.01 -0.30 -0.48
[5,] -0.21 -0.14  0.08  0.23  1.00  0.18  0.00 -0.26
[6,] -0.22 -0.35 -0.26 -0.01  0.18  1.00  0.20 -0.03
[7,]  0.11 -0.20 -0.42 -0.30  0.00  0.20  1.00  0.26
[8,]  0.28  0.08 -0.27 -0.48 -0.26 -0.03  0.26  1.00
> plot(factor.pa(circ.sim(16,500),2)$loadings,main="Circumplex Structure") #circumplex structure
> #
> #
> plot(factor.pa(item.sim(16,500),2)$loadings,main="Simple Structure") #simple structure
> #
> cluster.plot(factor.pa(item.dichot(16,low=0,high=1),2))
> 
> 
> 
> cleanEx(); nameEx("sim.structural")
> ### * sim.structural
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.structure
> ### Title: Create correlation matrices or data matrices with a particular
> ###   measurement and structural model
> ### Aliases: make.structural sim.structure sim.structural
> ### Keywords: multivariate datagen
> 
> ### ** Examples
> 
> fx <-matrix(c( .9,.8,.6,rep(0,4),.6,.8,-.7),ncol=2)              
> fy <- c(.6,.5,.4)
> Phi <-matrix( c(1,0,.7,.0,1,.7,.7,.7,1),ncol=3)
> gre.gpa <- sim.structural(fx,Phi,fy)
Loading required package: MASS
> print(gre.gpa,2)  
Call: NULL

 $model (Population correlation matrix) 
     V1   V2    V3    V4    V5    V6    V7    V8
V1 1.00 0.72  0.54  0.00  0.00  0.38  0.32  0.25
V2 0.72 1.00  0.48  0.00  0.00  0.34  0.28  0.22
V3 0.54 0.48  1.00  0.48 -0.42  0.50  0.42  0.34
V4 0.00 0.00  0.48  1.00 -0.56  0.34  0.28  0.22
V5 0.00 0.00 -0.42 -0.56  1.00 -0.29 -0.24 -0.20
V6 0.38 0.34  0.50  0.34 -0.29  1.00  0.30  0.24
V7 0.32 0.28  0.42  0.28 -0.24  0.30  1.00  0.20
V8 0.25 0.22  0.34  0.22 -0.20  0.24  0.20  1.00

$reliability (population reliability) 
 Vx1  Vx2  Vx3  Vx4  Vx5  Vy1  Vy2  Vy3 
0.81 0.64 0.72 0.64 0.49 0.36 0.25 0.16 
> round(correct.cor(gre.gpa$model,gre.gpa$reliability),2)  #correct for attenuation to see structure
     V1   V2    V3    V4    V5    V6    V7    V8
V1 0.81 1.00  0.71  0.00  0.00  0.70  0.70  0.70
V2 0.72 0.64  0.71  0.00  0.00  0.70  0.70  0.70
V3 0.54 0.48  0.72  0.71 -0.71  0.99  0.99  0.99
V4 0.00 0.00  0.48  0.64 -1.00  0.70  0.70  0.70
V5 0.00 0.00 -0.42 -0.56  0.49 -0.70 -0.70 -0.70
V6 0.38 0.34  0.50  0.34 -0.29  0.36  1.00  1.00
V7 0.32 0.28  0.42  0.28 -0.24  0.30  0.25  1.00
V8 0.25 0.22  0.34  0.22 -0.20  0.24  0.20  0.16
> 
> congeneric <- sim.structural(f=c(.9,.8,.7,.6)) # a congeneric model 
> congeneric 
Call: NULL

 $model (Population correlation matrix) 
     V1   V2   V3   V4
V1 1.00 0.72 0.63 0.54
V2 0.72 1.00 0.56 0.48
V3 0.63 0.56 1.00 0.42
V4 0.54 0.48 0.42 1.00

$reliability (population reliability) 
[1] 0.81 0.64 0.49 0.36
> 
> 
> 
> 
> cleanEx(); nameEx("simulation.circ")
> ### * simulation.circ
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: simulation.circ
> ### Title: Simulations of circumplex and simple structure
> ### Aliases: simulation.circ circ.simulation simulation.circ circ.sim.plot
> ### Keywords: multivariate datagen
> 
> ### ** Examples
> 
> 
> demo <- simulation.circ()
>  boxplot(demo[3:14])
>  title("4 tests of Circumplex Structure",sub="Circumplex, Ellipsoid, Simple Structure")
> circ.sim.plot(demo[3:14])  #compare these results to real data
> 
> 
> 
> 
> cleanEx(); nameEx("skew")
> ### * skew
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: skew
> ### Title: Calculate skew or kurtosis for a vector, matrix, or data.frame
> ### Aliases: skew kurtosi
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> round(skew(attitude),2)
[1] -0.36 -0.22  0.38 -0.05  0.20 -0.87  0.85
> round(kurtosi(attitude),2)
[1] -0.77 -0.68 -0.41 -1.22 -0.60  0.17  0.47
> 
> 
> 
> 
> cleanEx(); nameEx("smc")
> ### * smc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: smc
> ### Title: Find the Squared Multiple Correlation (SMC) of each variable
> ###   with the remaining variables in a matrix
> ### Aliases: smc SMC
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> R <- make.hierarchical()
Loading required package: MASS
> round(smc(R),2)
  V1   V2   V3   V4   V5   V6   V7   V8   V9 
0.44 0.37 0.28 0.30 0.24 0.17 0.18 0.14 0.09 
>  
> 
> 
> cleanEx(); nameEx("structure.graph")
> ### * structure.graph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: structure.graph
> ### Title: Draw a structural equation model specified by two measurement
> ###   models and a structural model
> ### Aliases: structure.graph structure.sem
> ### Keywords: multivariate hplot
> 
> ### ** Examples
> 
> fx <- matrix(c(.9,.8,.6,rep(0,4),.6,.8,-.7),ncol=2)
> fy <- matrix(c(.6,.5,.4),ncol=1)
> Phi <- matrix(c(1,0,0,0,1,0,.7,.7,1),ncol=3,byrow=TRUE)
> if(require(Rgraphviz)) { f1 <- structure.graph(fx,Phi,fy)  } else {f1 <- structure.sem(fx,Phi,fy)}
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
> 
> #symbolic input
> X2 <- matrix(c("a",0,0,"b","e1",0,0,"e2"),ncol=4)
> colnames(X2) <- c("X1","X2","E1","E2")
> phi2 <- diag(1,4,4)
> phi2[2,1] <- phi2[1,2] <- "r"
> if(require(Rgraphviz)) { f2 <- structure.graph(X2,Phi=phi2,errors=FALSE)  } else {f2 <- structure.sem(X2,Phi=phi2,errors=FALSE) }
> 
> #symbolic input with error 
> X2 <- matrix(c("a",0,0,"b"),ncol=2)
> colnames(X2) <- c("X1","X2")
> phi2 <- diag(1,2,2)
> phi2[2,1] <- phi2[1,2] <- "r"
> if(require(Rgraphviz)) { f3 <- structure.graph(X2,Phi=phi2)  } else {f3 <- structure.sem(X2,Phi=phi2)}
> 
> #and yet another one
> X6 <- matrix(c("a","b","c",rep(0,6),"d","e","f"),nrow=6)
> colnames(X6) <- c("L1","L2")
> rownames(X6) <- c("x1","x2","x3","x4","x5","x6")
> Y3 <- matrix(c("u","w","z"),ncol=1)
> colnames(Y3) <- "Y"
> rownames(Y3) <- c("y1","y2","y3")
> phi21 <- matrix(c(1,0,"r1",0,1,"r2",0,0,1),ncol=3)
> colnames(phi21) <- rownames(phi21) <-  c("L1","L2","Y")
> if(require(Rgraphviz)) {f4 <- structure.graph(X6,phi21,Y3,title="Structural model")} else {f4 <- structure.sem(X6,phi21,Y3)}
> 
> # and finally, a regression model
> X7 <- matrix(c("a","b","c","d","e","f"),nrow=6)
> if(require(Rgraphviz)) {f5 <- structure.graph(X7,regression=TRUE)}
> 
> #and a really messy regession model
> x8 <- c("b1","b2","b3")
> r8 <- matrix(c(1,"r12","r13","r12",1,"r23","r13","r23",1),ncol=3)
> if(require(Rgraphviz)) {f6<- structure.graph(x8,Phi=r8,regression=TRUE)}
> 
> 
> 
> cleanEx(); nameEx("structure.list")
> ### * structure.list
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: structure.list
> ### Title: Create factor model matrices from an input list
> ### Aliases: structure.list phi.list
> ### Keywords: multivariate models
> 
> ### ** Examples
> 
> fx <- structure.list(9,list(F1=c(1,2,3),F2=c(4,5,6),F3=c(7,8,9)))
> fy <- structure.list(3,list(Y=c(1,2,3)),"Y")
> phi <- phi.list(4,list(F1=c(4),F2=c(1,4),F3=c(2),F4=c(1,2,3)))
> fx
      F1   F2   F3  
 [1,] "a1" "0"  "0" 
 [2,] "a2" "0"  "0" 
 [3,] "a3" "0"  "0" 
 [4,] "0"  "b4" "0" 
 [5,] "0"  "b5" "0" 
 [6,] "0"  "b6" "0" 
 [7,] "0"  "0"  "c7"
 [8,] "0"  "0"  "c8"
 [9,] "0"  "0"  "c9"
> phi
   F1    F2    F3    F4   
F1 "1"   "rba" "0"   "rda"
F2 "0"   "1"   "rcb" "rdb"
F3 "0"   "0"   "1"   "rdc"
F4 "rad" "rbd" "0"   "1"  
> fy
     Y    
[1,] "Ya1"
[2,] "Ya2"
[3,] "Ya3"
> 
> 
> 
> 
> cleanEx(); nameEx("super.matrix")
> ### * super.matrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: super.matrix
> ### Title: Form a super matrix from two sub matrices.
> ### Aliases: super.matrix
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> mx <- matrix(c(.9,.8,.7,rep(0,4),.8,.7,.6),ncol=2)
> my <- matrix(c(.6,.5,.4))
> colnames(mx) <- paste("X",1:dim(mx)[2],sep="")
> rownames(mx) <- paste("Xv",1:dim(mx)[1],sep="")
> colnames(my) <- "Y"
> rownames(my) <- paste("Yv",1:3,sep="")
> super.matrix(mx,my)
     X1  X2   Y
Xv1 0.9 0.0 0.0
Xv2 0.8 0.0 0.0
Xv3 0.7 0.8 0.0
Xv4 0.0 0.7 0.0
Xv5 0.0 0.6 0.0
Yv1 0.0 0.0 0.6
Yv2 0.0 0.0 0.5
Yv3 0.0 0.0 0.4
> 
> 
> 
> cleanEx(); nameEx("table2matrix")
> ### * table2matrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: table2matrix
> ### Title: Convert a table with counts to a matrix or data.frame
> ###   representing those counts.
> ### Aliases: table2matrix table2df
> ### Keywords: models
> 
> ### ** Examples
> 
> data(cubits)
> cubit <- table2matrix(cubits,labs=c("height","cubit"))
> describe(cubit)
       var   n  mean   sd median trimmed  mad  min   max range  skew kurtosis
height   1 348 67.07 2.36  67.00   67.11 2.97 63.0 71.00  8.00 -0.09    -0.92
cubit    2 348 18.10 0.78  18.25   18.11 0.74 16.5 19.75  3.25 -0.09    -0.60
         se
height 0.13
cubit  0.04
> ellipses(cubit,n=1)
> 
> 
> 
> cleanEx(); nameEx("test.psych")
> ### * test.psych
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: test.psych
> ### Title: Testing of functions in the psych package
> ### Aliases: test.psych
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> test <- test.psych()
too many factors requested for this number of variables to use SMC, 1s used instead
Loading required package: GPArotation
In omega, 3 factors are too many for 4 variables using mle.  Using minres instead
In fa, too many factors requested for this number of variables to use SMC for communality estimates, 1s are used instead
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: grid
In fa, too many factors requested for this number of variables to use SMC for communality estimates, 1s are used instead
In smc, the correlation matrix was not invertible, smc's returned as 1s
Warning in log(det(m.inv.r)) : NaNs produced
In factor.stats, the correlation matrix is singular, an approximation is used
In factor.scores, the correlation matrix is singular, an approximation is used
Warning in log(det(m.inv.r)) : NaNs produced
In factor.stats, the correlation matrix is singular, an approximation is used
In factor.scores, the correlation matrix is singular, an approximation is used
Loading required package: MASS
> 
> 
> 
> 
> cleanEx(); nameEx("thurstone")
> ### * thurstone
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: thurstone
> ### Title: Thurstone Case V scaling
> ### Aliases: thurstone
> ### Keywords: models
> 
> ### ** Examples
> 
> data(vegetables)
> thurstone(veg)
Thurstonian scale (case 5) scale values 
Call: thurstone(x = veg)
   Turn     Cab    Beet     Asp     Car    Spin S.Beans    Peas    Corn 
   0.00    0.52    0.65    0.98    1.12    1.14    1.40    1.44    1.63 

 Goodness of fit of model   0.99> 
> 
> 
> cleanEx(); nameEx("tr")
> ### * tr
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: tr
> ### Title: Find the trace of a square matrix
> ### Aliases: tr
> ### Keywords: multivariate
> 
> ### ** Examples
> 
>  m <- matrix(1:16,ncol=4)
>  m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
>  tr(m)
[1] 34
>  
> 
> 
> 
> cleanEx(); nameEx("vegetables")
> ### * vegetables
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: vegetables
> ### Title: Paired comparison of preferences for 9 vegetables
> ### Aliases: vegetables veg
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(vegetables)
> thurstone(veg)
Thurstonian scale (case 5) scale values 
Call: thurstone(x = veg)
   Turn     Cab    Beet     Asp     Car    Spin S.Beans    Peas    Corn 
   0.00    0.52    0.65    0.98    1.12    1.14    1.40    1.44    1.63 

 Goodness of fit of model   0.99> 
> 
> 
> cleanEx(); nameEx("winsor")
> ### * winsor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: winsor
> ### Title: Find the Winsorized scores or means for a vector, matrix, or
> ###   data.frame
> ### Aliases: winsor winsor.means
> ### Keywords: univar
> 
> ### ** Examples
> 
> data(sat.act)
> winsor.means(sat.act) #compare with the means of the winsorized scores
    gender  education        age        ACT       SATV       SATQ 
  1.647143   3.391429  23.954286  28.957143 615.570000 614.521106 
> y <- winsor(sat.act)
> describe(y)
          var   n   mean    sd median trimmed    mad   min max range  skew
gender      1 700   1.65  0.48      2    1.68   0.00   1.0   2   1.0 -0.61
education   2 700   3.39  1.03      3    3.36   1.48   2.0   5   3.0  0.27
age         3 700  23.95  5.11     22   23.57   4.45  19.0  32  13.0  0.56
ACT         4 700  28.96  3.18     29   28.97   4.45  24.8  33   8.2 -0.06
SATV        5 700 615.57 72.79    620  618.21 118.61 510.0 700 190.0 -0.24
SATQ        6 687 614.52 80.88    620  616.87 118.61 500.0 710 210.0 -0.24
          kurtosis   se
gender       -1.62 0.02
education    -1.07 0.04
age          -1.30 0.19
ACT          -1.56 0.12
SATV         -1.43 2.75
SATQ         -1.47 3.09
> xy <- data.frame(sat.act,y)
> pairs.panels(xy) #to see the effect of winsorizing 
> x <- matrix(1:100,ncol=5)
> winsor(x)
      [,1] [,2] [,3] [,4] [,5]
 [1,]  4.8 24.8 44.8 64.8 84.8
 [2,]  4.8 24.8 44.8 64.8 84.8
 [3,]  4.8 24.8 44.8 64.8 84.8
 [4,]  4.8 24.8 44.8 64.8 84.8
 [5,]  5.0 25.0 45.0 65.0 85.0
 [6,]  6.0 26.0 46.0 66.0 86.0
 [7,]  7.0 27.0 47.0 67.0 87.0
 [8,]  8.0 28.0 48.0 68.0 88.0
 [9,]  9.0 29.0 49.0 69.0 89.0
[10,] 10.0 30.0 50.0 70.0 90.0
[11,] 11.0 31.0 51.0 71.0 91.0
[12,] 12.0 32.0 52.0 72.0 92.0
[13,] 13.0 33.0 53.0 73.0 93.0
[14,] 14.0 34.0 54.0 74.0 94.0
[15,] 15.0 35.0 55.0 75.0 95.0
[16,] 16.0 36.0 56.0 76.0 96.0
[17,] 16.2 36.2 56.2 76.2 96.2
[18,] 16.2 36.2 56.2 76.2 96.2
[19,] 16.2 36.2 56.2 76.2 96.2
[20,] 16.2 36.2 56.2 76.2 96.2
> winsor.means(x)
[1] 10.5 30.5 50.5 70.5 90.5
> y <- 1:11
> winsor(y,trim=.5)
 [1] 6 6 6 6 6 6 6 6 6 6 6
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed:  66.595 1.041 69.245 0.003 0.003 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')