tests/Thurstone.R

#Example from: Gradient Projection Algorithms and Software for 
#  Arbitrary Rotation Criteria in Factor Analysis.
#  by Coen A. Bernaards and Robert I. Jennrich
#  Website: http://www.stat.ucla.edu/research

 Sys.getenv("R_LIBS")
 library()
 require("GPArotation")
 search()
 Sys.info()

 data("Thurstone", package="GPArotation") 
 if (!exists("box20")) stop("Test data not found. Testing stopped.")

fuzz <- 1e-5 
all.ok <- TRUE  

#  Thurstone's box problem.  (1947, p. 136) 
# The matrix box20 is the initial loading matrix from Thurstone's box problem.
            
# This takes a lot of iterations to converge at a higher tolerance
qbox20 <- quartimax(box20, eps=1e-5)
qbox20G <- GPForth(box20, Tmat=diag(1,3), method="quartimax", eps=1e-5)

 if( fuzz < max(abs(qbox20$loadings - qbox20G$loadings))) {
    cat("Calculated value is not the same as test value in test Thurstone 1. Value:\n")
    print(qbox20$loadings - qbox20G$loadings, digits=18)
    cat("difference:\n")
    print(qbox20$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 

#qbox20$Th - qbox20G$Th 

# These values compare with those in:
#    http://www.stat.ucla.edu/research/web.pdf
 tst <- t(matrix(c(
   0.0104916072210123716, -0.993396087928394733, -0.089861775335686706,
   0.1584646383898045685, -0.167305085570175344, -0.967087879524061056,
   0.9822741057703969769, -0.094961339079248266, -0.081938545344928893,
   0.1249962020162782989, -0.597065497283680413, -0.789290657131387352,
   0.8695614167874907707, -0.471622450093366785, -0.090438968384549553,
   0.8757114893176747294, -0.141012080768234127, -0.452333925937943637,
   0.0679423211019681700, -0.811411071716238719, -0.588554936857709099,
   0.4066768108416509708, -0.907862149146695163, -0.115673202040957226,
   0.5770808894249742638, -0.142370726163931066, -0.806527261406603468,
   0.1012712863762783577, -0.723336747696182614, -0.694640249329106285,
   0.5000928657774492692, -0.949746569049947253, -0.046846346456817907,
   0.7412589798326677526, -0.140350561965914555, -0.663578062154924320,
   0.0055655501003109590, -0.983847100401775698, -0.120037109608235590,
   0.2142330103903098415, -0.119429100752156334, -0.947421187831809397,
   0.9550804066106526324, -0.108275659756619305, -0.039227521113362487,
   0.7823218737697450464, -0.405437596810190704, -0.439275358874331168,
   0.3626971102221024923, -0.753122462957226402, -0.546281394544768872,
   0.0162483298780003657, -0.966230359337758582, -0.052114148464710915,
   0.1076692386876715729, -0.206734953950642314, -0.934620775424686911,
   0.9744239420161749932, -0.092650552854598708, -0.090828719474599584
   ), 3, 20))
 
 if( fuzz < max(abs(qbox20$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 2. Value:\n")
    print(qbox20$loadings, digits=18)
    cat("difference:\n")
    print(qbox20$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 

 tst <- t(matrix(c(
   0.57232345894276127, -0.60751194947821441, -0.55079496147384377,
   0.60249460283341838,  0.76716797198365361, -0.22012168525406509,
   0.55627880770383020, -0.20587018726291534,  0.80509089803322043
   ), 3, 3))
 
 if( fuzz < max(abs(qbox20$Th - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 3. Value:\n")
    print(qbox20$Th, digits=18)
    cat("difference:\n")
    print(qbox20$Th - tst, digits=18)
    all.ok <- FALSE  
    } 


# sorted absolute loading plots. 
sal <- abs(c(loadings(qbox20)))[order(abs(c(loadings(qbox20))))]
plot(seq(length(sal)), sal)


 #compare quartimax rotation of the initial loading matrix box20.

 if( fuzz < max(abs(loadings(qbox20) - box20 %*% qbox20$Th ))) {
    cat("Calculated value is not the same as test value in test Thurstone 4. Value:\n")
    print(loadings(qbox20), digits=18)
    cat("difference:\n")
    print(loadings(qbox20) - box20 %*% qbox20$Th, digits=18)
    all.ok <- FALSE  
    } 


            
qminbox20G <- GPFoblq(box20, Tmat=diag(1,3), method="quartimin", eps=1e-5)
qminbox20  <- quartimin(box20, eps=1e-5)

 if( fuzz < max(abs(loadings(qminbox20) - qminbox20G$loadings))) {
    cat("Calculated value is not the same as test value in test Thurstone 5. Value:\n")
    print(qminbox20G$loadings , digits=18)
    cat("difference:\n")
    print(loadings(qminbox20) - qminbox20G$loadings, digits=18)
    all.ok <- FALSE  
    } 

#qminbox20$Th - quartimin(box20)$Th 

# These values compare with those in:
#    http://www.stat.ucla.edu/research/web.pdf
 tst <- t(matrix(c(
   -0.099561899210599963, -1.0236437309424475384,  0.017110338313848200,
   -0.007103778102200991,  0.0427848301281630802, -1.009962780073245581,
    1.012864497258948226,  0.0331727792925069487,  0.050367710973030555,
   -0.054843850612513692, -0.4493155290974688021, -0.772334543778026350,
    0.856287122381722998, -0.3740197232441037078,  0.069350368268248391,
    0.835580575619599641,  0.0487450425576793633, -0.360381644212301344,
   -0.102893671454670210, -0.7226715938020771279, -0.537456650126404090,
    0.322103633211960838, -0.8816846447967544576,  0.031159743715387874,
    0.462799683447739529,  0.0852338438217692257, -0.783762970578423479,
   -0.076585435689138226, -0.6043060025891554554, -0.658295846696152820,
    0.427772530893690217, -0.9288687512327726825,  0.122866182561916254,
    0.659408232467282085,  0.0772080094990600374, -0.607348040513722709,
   -0.108761719100651882, -1.0079608432113262850, -0.017378089000366713,
    0.059518597564186392,  0.0955950614351480238, -0.986779686330629513,
    0.989890866913205381,  0.0071520817823045348,  0.094691644950703049,
    0.713733277219835149, -0.2427293600063723522, -0.328268187306521242,
    0.220344503737931546, -0.6353746612195683152, -0.459661643730432223,
   -0.084703580704062989, -1.0022284232457450148,  0.055740317456252478,
   -0.059151779416785115, -0.0113377397453605679, -0.976867596293413132,
    1.003360458549731771,  0.0365098037316876067,  0.039427150580815938
    ), 3, 20))
 
 if( fuzz < max(abs(qminbox20G$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 6. Value:\n")
    print(qminbox20G$loadings, digits=18)
    cat("difference:\n")
    print(qminbox20G$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 

 tst <- t(matrix(c(
   1.00000000000000000, -0.25676300454795098, -0.32155119431295237,
  -0.25676300454795098,  1.00000000000000000,  0.33656790396842257,
  -0.32155119431295237,  0.33656790396842257,  1.00000000000000000
   ), 3, 3))
 
 if( fuzz < max(abs(qminbox20G$Phi - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 7. Value:\n")
    print(qminbox20G$Phi, digits=18)
    cat("difference:\n")
    print(qminbox20G$Phi - tst, digits=18)
    all.ok <- FALSE  
    } 


#To fuzz precision the rotated loading matrix and the factor cor-
#relation matrix phi are identical to those produced using the oblique GP
#algorithm with exact derivatives.

 if( fuzz < max(abs(qminbox20G$Phi -
                  t(qminbox20G$Th )%*% qminbox20G$Th ))) {
    cat("Calculated value is not the same as test value in test Thurstone 8. Value:\n")
    print(qminbox20G$Phi, digits=18)
    cat("difference:\n")
    print(qminbox20G$Phi - 
                  t(qminbox20G$Th )%*% qminbox20G$Th, digits=18)
    all.ok <- FALSE  
    } 

 #compare quartimin rotation of the initial loading matrix box20.
   if( fuzz < max(abs(qminbox20G$loadings - box20 %*% solve(t(qminbox20G$Th))))) {
    cat("Calculated value is not the same as test value in test Thurstone 9. Value:\n")
    print(qminbox20G$loadings, digits=18)
    cat("difference:\n")
    print(qminbox20G$loadings - box20 %*% solve(t(qminbox20G$Th)), digits=18)
    all.ok <- FALSE  
    } 

 data("box26", package="GPArotation") 
 if (!exists("box26")) stop("Test data box26 not found. Testing stopped.")

qbox26  <- GPForth(box26, Tmat=diag(1,3), method="quartimax", eps=1e-5)

 tst <- t(matrix(c(
   0.6245197355925140581, -0.2708954695931116152,  0.7151983951389878635,
   0.7386116884036847408,  0.6266342260884526505, -0.0617439911892987553,
   0.7803093788467402314, -0.3830982859243221017, -0.4578886072022986253,
   0.8540550453155928423,  0.2886436985992582027,  0.4062915145925659610,
   0.8810593765418006651, -0.4428658074662961130,  0.1233946983666596581,
   0.9084731768740617053,  0.1540526132602804965, -0.3723026715563940159,
   0.8150592858039771293,  0.0441965358534676597,  0.5600768044145943980,
   0.8466584455973064083,  0.4551177395514792168,  0.1889929089788950356,
   0.8156808837280125069, -0.4090629943132625956,  0.3690652552112651530,
   0.9629492340906220527, -0.4781483041690369196, -0.0866081507974762743,
   0.8731366884896356595,  0.3451069860590937899, -0.2914969834947889749,
   0.8921854600753849063, -0.0276323108621970258, -0.4257376659710629951,
  -0.0938760381595044741, -0.7873218033841372643,  0.6012450975895150540,
   0.0938760381595044741,  0.7873218033841372643, -0.6012450975895150540,
  -0.0986092863860908303,  0.1513605567468480073,  0.9692559984337008050,
   0.0986092863860908303, -0.1513605567468480073, -0.9692559984337008050,
  -0.0189573629854957251,  0.9527983290277913797,  0.2944078167958268377,
   0.0189573629854957251, -0.9527983290277913797, -0.2944078167958268377,
   0.8394181189595459891,  0.3631767908642606346,  0.3398717995655929913,
   0.8703065201362156778, -0.4691145408161159214,  0.0770980453920554615,
   0.9141063746617547059,  0.1583184861345137973, -0.3535658252020681958,
   0.8348118627305495254,  0.3535663452183119837,  0.3271666140872500073,
   0.8541352373790773722, -0.4476738735312740247,  0.0569042988261160704,
   0.9034738474019414767,  0.1663655738425987851, -0.3227406124130587362,
   0.9861758757457432800,  0.0103496363116840455,  0.0635926656567585569,
   0.9643516568468981642,  0.0660181478622221818, -0.0304218028637989850
   ), 3, 26))
 
 if( fuzz < max(abs(qbox26$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 10. Value:\n")
    print(qbox26$loadings, digits=18)
    cat("difference:\n")
    print(qbox26$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 



 tst <- t(matrix(c(
   0.9996572020207266096, 0.0216275672176080257,  0.0147555679097727491,
  -0.0158190757965277796, 0.9480178905874908635, -0.3178235925273457108,
  -0.0208622934749700742, 0.3174812237948764770,  0.9480350400953921897
   ), 3, 3))
 
 if( fuzz < max(abs(qbox26$Th - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 11. Value:\n")
    print(qbox26$Th, digits=18)
    cat("difference:\n")
    print(qbox26$Th - tst, digits=18)
    all.ok <- FALSE  
    } 


qminbox26 <- GPFoblq(box26, Tmat=diag(1,3), method="quartimin", eps=1e-5)

 tst <- t(matrix(c(
   0.6088436426802223966, -0.2567107018725688361,  0.7213648290819488773,
   0.7318447535507376367,  0.6298398026581654152, -0.0549983771960348838,
   0.7973321695017724364, -0.3855960314746548212, -0.4504478973568259437,
   0.8392144987741166906,  0.2994932968625432235,  0.4143558581243267924,
   0.8833452352200144020, -0.4361046712803113290,  0.1319331147095905710,
   0.9161366872228343672,  0.1535557844336666034, -0.3638337328539109072,
   0.7993355454002614158,  0.0571270784641514026,  0.5678963531379384033,
   0.8354288250614068101,  0.4626764152757318893,  0.1968789749765105790,
   0.8109923806202916641, -0.3989909333845649830,  0.3770226870580207779,
   0.9712737747877250305, -0.4740722765307348041, -0.0773243882106463137,
   0.8761501947960563808,  0.3456235893514668089, -0.2834183138879167174,
   0.9036601763684347643, -0.0290211959776035672, -0.4173652812159966974,
  -0.0995525797764766768, -0.7788574612781464790,  0.6007791331268093060,
   0.0995525797764766768,  0.7788574612781464790, -0.6007791331268093060,
  -0.1264036712449473909,  0.1653130238928011975,  0.9684661160120416890,
   0.1264036712449473909, -0.1653130238928011975, -0.9684661160120416890,
  -0.0392946742598458687,  0.9571059478962877787,  0.2939285303852590125,
   0.0392946742598458687, -0.9571059478962877787, -0.2939285303852590125,
   0.8253744379910458173,  0.3729516010405902748,  0.3477554718030251846,
   0.8741734789142978634, -0.4631063486451737488,  0.0855349365396926159,
   0.9212130243051569467,  0.1581334796580046720, -0.3450412531516501846,
   0.8212340853954427367,  0.3631252613622908965,  0.3350076577679809153,
   0.8582635618771776720, -0.4420579024138228674,  0.0651757040961165046,
   0.9096561314838297330,  0.1665824736284239604, -0.3143133889989875307,
   0.9840845767693481294,  0.0168070160966761091,  0.0729425956763933708,
   0.9640420478016114014,  0.0709475796833391181, -0.0213192081807395371
   ), 3, 26))
 
 if( fuzz < max(abs(qminbox26$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 12. Value:\n")
    print(qminbox26$loadings, digits=18)
    cat("difference:\n")
    print(qminbox26$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 



 tst <- t(matrix(c(
   1.000000000000000	 ,  0.00767934084449363279,  0.0170654511973979163,
   0.00767934084449363279,  1.000000000000000	  , -0.0144994900961642244,
   0.01706545119739791630, -0.01449949009616422445,  1.000000000000000     
    ), 3, 3))
 
 if( fuzz < max(abs(qminbox26$Phi - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 13. Value:\n")
    print(qminbox26$Phi, digits=18)
    cat("difference:\n")
    print(qminbox26$Phi - tst, digits=18)
    all.ok <- FALSE  
    } 



 tst <- t(matrix(c(
    0.9993401424148040668, 0.0347479564402226465,  0.0408645923859655008,
   -0.0179660947915933414, 0.9476477730670300748, -0.3324117322929439067,
   -0.0315673755054017430, 0.3174212937474846785,  0.9422486536594960604
    ), 3, 3))
 
 if( fuzz < max(abs(qminbox26$Th - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 14. Value:\n")
    print(qminbox26$Th, digits=18)
    cat("difference:\n")
    print(qminbox26$Th - tst, digits=18)
    all.ok <- FALSE  
    } 



cat("tests completed.\n")

if (! all.ok) stop("some tests FAILED")


cat("tests completed.\n")

if (! all.ok) stop("some tests FAILED")

Try the GPArotation package in your browser

Any scripts or data that you put into this service are public.

GPArotation documentation built on April 16, 2022, 5:05 p.m.