R Under development (unstable) (2024-02-02 r85855 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
>  Sys.getenv("R_LIBS")
[1] "D:\\temp\\RtmpO4cY16\\RLIBS_2b7c074005b36"
>  library()
>  require("GPArotation")
Loading required package: GPArotation
>  search()
 [1] ".GlobalEnv"          "package:GPArotation" "package:stats"      
 [4] "package:graphics"    "package:grDevices"   "package:utils"      
 [7] "package:datasets"    "package:methods"     "Autoloads"          
[10] "package:base"       
>  Sys.info()
       sysname        release        version       nodename        machine 
     "Windows"   "Server x64"  "build 20348"     "CRANWIN3"       "x86-64" 
         login           user effective_user 
        "CRAN"         "CRAN"         "CRAN" 
> 
> #require("stats")  
> 
> fuzz <- 1e-6 
> all.ok <- TRUE  
> 
> 
> # test MASS 4th ed. p 322-324
> 
>   data(ability.cov)
>   ability.cov
$cov
        general picture  blocks   maze reading   vocab
general  24.641   5.991  33.520  6.023  20.755  29.701
picture   5.991   6.700  18.137  1.782   4.936   7.204
blocks   33.520  18.137 149.831 19.424  31.430  50.753
maze      6.023   1.782  19.424 12.711   4.757   9.075
reading  20.755   4.936  31.430  4.757  52.604  66.762
vocab    29.701   7.204  50.753  9.075  66.762 135.292

$center
[1] 0 0 0 0 0 0

$n.obs
[1] 112

>   ability.FA <- factanal(factors = 1, covmat=ability.cov)
> 
>   (ability.FA <- update(ability.FA, factors = 2))

Call:
factanal(factors = 2, covmat = ability.cov)

Uniquenesses:
general picture  blocks    maze reading   vocab 
  0.455   0.589   0.218   0.769   0.052   0.334 

Loadings:
        Factor1 Factor2
general 0.499   0.543  
picture 0.156   0.622  
blocks  0.206   0.860  
maze    0.109   0.468  
reading 0.956   0.182  
vocab   0.785   0.225  

               Factor1 Factor2
SS loadings      1.858   1.724
Proportion Var   0.310   0.287
Cumulative Var   0.310   0.597

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 6.11 on 4 degrees of freedom.
The p-value is 0.191 
> 
> #  ability.FA2 <- factanal(factors = 2, covmat = ability.cov)
> # max(abs(ability.FA2$loadings - ability.FA$loadings))
> 
> #  summary(ability.FA) MASS ed.4 p 323 seems to be print not summary in R 2.0.1
>   ability.FA

Call:
factanal(factors = 2, covmat = ability.cov)

Uniquenesses:
general picture  blocks    maze reading   vocab 
  0.455   0.589   0.218   0.769   0.052   0.334 

Loadings:
        Factor1 Factor2
general 0.499   0.543  
picture 0.156   0.622  
blocks  0.206   0.860  
maze    0.109   0.468  
reading 0.956   0.182  
vocab   0.785   0.225  

               Factor1 Factor2
SS loadings      1.858   1.724
Proportion Var   0.310   0.287
Cumulative Var   0.310   0.597

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 6.11 on 4 degrees of freedom.
The p-value is 0.191 
> 
> # this is default varimax rotation. There are 3rd+ digit differences with MASS
>   tst <- t(matrix(c(
+       0.499437829039896530, 0.54344904693111962,
+       0.156070079431279873, 0.62153798991197484,
+       0.205786989958578748, 0.85992588538426895,
+       0.108530754440558652, 0.46776101732283504,
+       0.956242470279811574, 0.18209631992182243,
+       0.784768183877880943, 0.22482213687364205
+       ), 2, 6))
>  
> 
>  if( fuzz < max(abs(loadings(ability.FA) - tst))) {
+     cat("Calculated value is not the same as test value in test 1. Value:\n")
+     #print(loadings(ability.FA), digits=18) this truncates
+     print(unclass(ability.FA$loadings), digits=18)
+     cat("difference:\n")
+     print(unclass(ability.FA$loadings) - tst, digits=18)
+     all.ok <- FALSE  
+     } 
> 
> 
>   # differences with MASS here are a bit more than might be expected,
>   # but there is already a difference before rotation.
>   (oblirot <- oblimin(loadings(ability.FA)))
Oblique rotation method Oblimin Quartimin converged.
Loadings:
        Factor1 Factor2
general  0.3864  0.4745
picture -0.0110  0.6459
blocks  -0.0263  0.8961
maze    -0.0180  0.4883
reading  0.9901 -0.0371
vocab    0.7906  0.0526

               Factor1 Factor2
SS loadings      1.825   1.757
Proportion Var   0.304   0.293
Cumulative Var   0.304   0.597

Phi:
        Factor1 Factor2
Factor1   1.000   0.465
Factor2   0.465   1.000
> 
>   obli2 <- factanal(factors = 2, covmat = ability.cov, rotation="oblimin")
> 
>   max(abs(loadings(oblirot) - loadings(obli2)))
[1] 2.38643e-06
> 
> 
> # factanal(factors = 2, covmat = ability.cov, scores = Bartlett, rotation="oblimin")
> 
> 
>  tst <- t(matrix(c(
+        0.3863637969729337152,  0.4745113977203344047,
+       -0.0110032278171669998,  0.6458708261423832253,
+       -0.0262888675561207576,  0.8961123879025085781,
+       -0.0180180060207963122,  0.4882918937716873575,
+        0.9900948712271664398, -0.0370729040114848238,
+        0.7905663749272058283,  0.0526099352008769991
+     ), 2, 6))
>  
>  if( fuzz < max(abs(loadings(oblirot) - tst ))) {
+     cat("Calculated value is not the same as test value in test 2. Value:\n")
+     print(loadings(oblirot), digits=18)
+     cat("difference:\n")
+     print(loadings(oblirot) - tst, digits=18)
+     all.ok <- FALSE  
+     } 
> 
> cat("tests completed.\n")
tests completed.
> 
> if (! all.ok) stop("some tests FAILED")
> 
> proc.time()
   user  system elapsed 
   0.17    0.07    0.25