pkgname <- "psych" source(file.path(R.home("share"), "R", "examples-header.R")) options(warn = 1) library('psych') base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') base::assign(".old_wd", base::getwd(), 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 #to test most of the psych package run the following #test.psych() cleanEx() nameEx("AUC") ### * AUC flush(stderr()); flush(stdout()) ### Name: AUC ### Title: Decision Theory measures of specificity, sensitivity, and d ### prime ### Aliases: AUC auc Specificity Sensitivity ### Keywords: multivariate ### ** Examples AUC(c(30,20,20,30)) #specify the table input AUC(c(140,60,100,900)) #Metz example with colors AUC(c(140,60,100,900),col=c("grey","black")) #Metz example 1 no colors AUC(c(80,120,40, 960)) #Metz example 2 Note how the accuracies are the same but d's differ AUC(c(49,40,79,336)) #Wiggins p 249 AUC(BR=.05,SR=.254,Phi = .317) #Wiggins 251 extreme Base Rates cleanEx() nameEx("Garcia") ### * Garcia flush(stderr()); flush(stdout()) ### Name: Garcia ### Title: Data from the sexism (protest) study of Garcia, Schmitt, ### Branscome, and Ellemers (2010) ### Aliases: Garcia protest GSBE ### Keywords: datasets ### ** Examples data(GSBE) #alias to Garcia data set ## Just do regressions with interactions setCor(respappr ~ prot2 * sexism,std=FALSE,data=Garcia,main="Moderated (mean centered )") setCor(respappr ~ prot2 * sexism,std=FALSE,data=Garcia,main="Moderated (don't center)", zero=FALSE) #demonstrate interaction plots plot(respappr ~ sexism, pch = 23- protest, bg = c("black","red", "blue")[protest], data=Garcia, main = "Response to sexism varies as type of protest") by(Garcia,Garcia$protest, function(x) abline(lm(respappr ~ sexism, data =x),lty=c("solid","dashed","dotted")[x$protest+1])) text(6.5,3.5,"No protest") text(3,3.9,"Individual") text(3,5.2,"Collective") #compare two models (bootstrapping n.iter set to 50 for speed # 1) mean center the variables prior to taking product terms mod1 <- mediate(respappr ~ prot2 * sexism +(sexism),data=Garcia,n.iter=50 ,main="Moderated mediation (mean centered)") # 2) do not mean center mod2 <- mediate(respappr ~ prot2 * sexism +(sexism),data=Garcia,zero=FALSE, n.iter=50, main="Moderated mediation (not centered") summary(mod1) summary(mod2) cleanEx() nameEx("Gleser") ### * Gleser flush(stderr()); flush(stdout()) ### Name: Gleser ### Title: Example data from Gleser, Cronbach and Rajaratnam (1965) to show ### basic principles of generalizability theory. ### Aliases: Gleser ### Keywords: datasets ### ** Examples #Find the MS for each component: #First, stack the data data(Gleser) stack.g <- stack(Gleser) st.gc.df <- data.frame(stack.g,Persons=rep(letters[1:12],12), Items=rep(letters[1:6],each=24),Judges=rep(letters[1:2],each=12)) #now do the ANOVA anov <- aov(values ~ (Persons*Judges*Items),data=st.gc.df) summary(anov) cleanEx() nameEx("Gorsuch") ### * Gorsuch flush(stderr()); flush(stdout()) ### Name: Gorsuch ### Title: Example data set from Gorsuch (1997) for an example factor ### extension. ### Aliases: Gorsuch ### Keywords: datasets ### ** Examples data(Gorsuch) Ro <- Gorsuch[1:6,1:6] Roe <- Gorsuch[1:6,7:10] fo <- fa(Ro,2,rotate="none") fa.extension(Roe,fo,correct=FALSE) cleanEx() nameEx("Harman") ### * Harman flush(stderr()); flush(stdout()) ### Name: Harman ### Title: Five data sets from Harman (1967). 9 cognitive variables from ### Holzinger and 8 emotional variables from Burt ### Aliases: Harman Harman.Burt Harman.Holzinger Harman.political Harman.5 ### Harman.8 ### Keywords: datasets ### ** Examples data(Harman) cor.plot(Harman.Holzinger) cor.plot(Harman.Burt) smc(Harman.Burt) #note how this produces impossible results f2 <- fa(Harman.8,2, rotate="none") #minres matches Harman and Jones 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) ICC(sf,lmer=FALSE) #just use the aov procedure #data(sai) sai <- psychTools::sai sai.xray <- subset(sai,(sai$study=="XRAY") & (sai$time==1)) xray.icc <- ICC(sai.xray[-c(1:3)],lmer=TRUE,check.keys=TRUE) xray.icc xray.icc$lme #show the variance components as well 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,title="ICLUST of the Harman data") summary(ic.out) #use all defaults and stop at 4 clusters ic.out4 <- iclust(test.data,nclusters =4,title="Force 4 clusters") summary(ic.out4) ic.out1 <- iclust(test.data,beta=3,beta.size=3) #use more stringent criteria ic.out #more complete output plot(ic.out4) #this shows the spatial representation #use a dot graphics viewer on the out.file #dot.graph <- ICLUST.graph(ic.out,out.file="test.ICLUST.graph.dot") #show the equivalent of a factor solution fa.diagram(ic.out4$pattern,Phi=ic.out4$Phi,main="Pattern taken from iclust") cleanEx() nameEx("ICLUST.graph") ### * ICLUST.graph flush(stderr()); flush(stdout()) ### Name: ICLUST.graph ### Title: create control code for ICLUST graphical output ### Aliases: ICLUST.graph 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{https://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 ic.out <- ICLUST(test.data) #uses iclust.diagram instead cleanEx() nameEx("KMO") ### * KMO flush(stderr()); flush(stdout()) ### Name: KMO ### Title: Find the Kaiser, Meyer, Olkin Measure of Sampling Adequacy ### Aliases: KMO ### Keywords: multivariate models ### ** Examples KMO(Thurstone) KMO(Harman.political) #compare to the results in Dziuban and Shirkey (1974) cleanEx() nameEx("Pinv") ### * Pinv flush(stderr()); flush(stdout()) ### Name: Pinv ### Title: Compute the Moore-Penrose Pseudo Inverse of a matrix ### Aliases: Pinv ### Keywords: multivariate ### ** Examples round(Pinv(Thurstone) %*% Thurstone,2) #an identity matrix sl <- schmid(Thurstone,3) #The schmid-leiman solution is less than full rank F <- sl$sl[,1:4] #the SL solution is general + 3 gropus R <- Thurstone # diag(R) <- sl$sl[,5] #the reproduced matrix (R - U2) S <- t(Pinv(t(F) %*% F) %*% t(F) %*% R) #the structure matrix Phi <- t(S) %*% F %*% Pinv(t(F) %*% F) #the factor covariances cleanEx() nameEx("Promax") ### * Promax flush(stderr()); flush(stdout()) ### Name: Promax ### Title: Perform Procustes,bifactor, promax or targeted rotations and ### return the inter factor angles. ### Aliases: Promax faRotate Procrustes TargetQ TargetT target.rot bifactor ### biquartimin varimin vgQ.bimin vgQ.targetQ vgQ.varimin equamax ### Keywords: multivariate models ### ** Examples jen <- sim.hierarchical() f3 <- fa(jen,3,rotate="varimax") f3 #not a very clean solution Promax(f3) #this obliquely rotates, but from the varimax target target.rot(f3) #this obliquely rotates to wards a simple structure target #compare this rotation with the solution from a targeted rotation aimed for #an independent cluster solution #now try a bifactor solution fb <-fa(jen,3,rotate="bifactor") fq <- fa(jen,3,rotate="biquartimin") #Suitbert Ertel has suggested varimin fm <- fa(jen,3,rotate="varimin") #the Ertel varimin fn <- fa(jen,3,rotate="none") #just the unrotated factors #compare them factor.congruence(list(f3,fb,fq,fm,fn)) # compare an oblimin with a target rotation using the Browne algorithm #note that we are changing the factor #order (this is for demonstration only) Targ <- make.keys(9,list(f1=1:3,f2=7:9,f3=4:6)) Targ <- scrub(Targ,isvalue=1) #fix the 0s, allow the NAs to be estimated Targ <- list(Targ) #input must be a list #show the target Targ fa(Thurstone,3,rotate="TargetQ",Target=Targ) #targeted oblique rotation #compare with oblimin f3 <- fa(Thurstone,3) #now try a targeted orthogonal rotation Targ <- make.keys(9,list(f1=1:3,f2=7:9,f3=4:6)) faRotate(f3$loadings,rotate="TargetT",Target=list(Targ)) #orthogonal 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) apply(attitude,2,sd) #all complete attitude[,1] <- NA SD(attitude) #missing a column describe(attitude) cleanEx() nameEx("Schmid.Leiman") ### * Schmid.Leiman flush(stderr()); flush(stdout()) ### Name: Schmid ### Title: 12 variables created by Schmid and Leiman to show the ### Schmid-Leiman Transformation ### Aliases: Schmid schmid.leiman West Chen ### Keywords: datasets ### ** Examples data(Schmid) cor.plot(Schmid,TRUE) print(fa(Schmid,6,rotate="oblimin"),cut=0) #shows an oblique solution round(cov2cor(schmid.leiman),2) cor.plot(cov2cor(West),TRUE) cleanEx() nameEx("Tucker") ### * Tucker flush(stderr()); flush(stdout()) ### Name: Tucker ### Title: 9 Cognitive variables discussed by Tucker and Lewis (1973) ### Aliases: Tucker ### Keywords: datasets ### ** Examples data(Tucker) fa(Tucker,2,n.obs=710) omega(Tucker,2) cleanEx() nameEx("VSS") ### * VSS flush(stderr()); flush(stdout()) ### Name: VSS ### Title: Apply the Very Simple Structure, MAP, and other criteria to ### determine the appropriate number of factors. ### Aliases: vss VSS MAP nfactors ### Keywords: multivariate models ### ** Examples #test.data <- Harman74.cor$cov #my.vss <- VSS(test.data,title="VSS of 24 mental tests") #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="minres" ,title="VSS of 24 circumplex variables") VSS(sim.item(nvar=24),fm="minres" ,title="VSS of 24 simple structure variables") 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 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 the successive eigen values for a scree test ### Aliases: VSS.scree scree ### Keywords: multivariate hplot ### ** Examples 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 tetrachoric, recreate table the table ### to create the Yule coefficient. ### Aliases: Yule Yule.inv Yule2phi Yule2tetra Yule2poly YuleBonett YuleCor ### Keywords: multivariate models ### ** Examples Nach <- matrix(c(40,10,20,50),ncol=2,byrow=TRUE) Yule(Nach) Yule.inv(.81818,c(50,60),n=120) Yule2phi(.81818,c(50,60),n=120) Yule2tetra(.81818,c(50,60),n=120) phi(Nach) #much less #or express as percents and do not specify n Nach <- matrix(c(40,10,20,50),ncol=2,byrow=TRUE) Nach/120 Yule(Nach) Yule.inv(.81818,c(.41667,.5)) Yule2phi(.81818,c(.41667,.5)) Yule2tetra(.81818,c(.41667,.5)) phi(Nach) #much less YuleCor(psychTools::ability[,1:4],,TRUE) YuleBonett(Nach,1) #Yule Q YuleBonett(Nach,.5) #Yule Y YuleBonett(Nach,.75) #Digby H YuleBonett(Nach,,TRUE) #Yule* is a generalized Yule 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 alpha.ci ### Keywords: models multivariate ### ** Examples set.seed(42) #keep the same starting values #four congeneric measures r4 <- sim.congeneric() alpha(r4) #nine hierarchical measures -- should actually use omega r9 <- sim.hierarchical() alpha(r9) # examples 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) #specify which items to reverse key by name alpha(two.f,keys=c("V3","V4","V5","V6")) cov.two <- cov(two.f) alpha(cov.two,check.keys=TRUE) #automatic reversal base upon first component alpha(two.f,check.keys=TRUE) #note that the median is much less than the average R #this suggests (correctly) that the 1 factor model is probably wrong #an example with discrete item responses -- show the frequencies items <- sim.congeneric(N=500,short=FALSE,low=-2,high=2, categorical=TRUE) #500 responses to 4 discrete items with 5 categories a4 <- alpha(items$observed) #item response analysis of congeneric measures a4 #summary just gives Alpha summary(a4) cleanEx() nameEx("anova.psych") ### * anova.psych flush(stderr()); flush(stdout()) ### Name: anova.psych ### Title: Model comparison for regression, mediation, and factor analysis ### Aliases: anova.psych ### Keywords: models multivariate ### ** Examples if(require("psychTools")) { m1 <- setCor(reaction ~ import, data = Tal_Or,std=FALSE) m2 <- setCor(reaction ~ import+pmi, data = Tal_Or,std=FALSE) m3 <- setCor(reaction ~ import+pmi + cond, data = Tal_Or,std=FALSE) anova(m1,m2,m3) } #Several interesting test cases are taken from analyses of the Spengler data set #Although the sample sizes are actually very large in the first wave, I use the #sample sizes from the last wave #This data set is actually in psychTools but is copied here until we can update psychTools #We set the n.iter to be 50 instead of the default value of 5,000 if(require("psychTools")) { mod1 <- mediate(Income.50 ~ IQ + Parental+ (Ed.11) ,data=Spengler, n.obs = 1952, n.iter=50) mod2 <- mediate(Income.50 ~ IQ + Parental+ (Ed.11) + (Income.11) ,data=Spengler,n.obs = 1952, n.iter=50) #Now, compare these models anova(mod1,mod2) } f3 <- fa(Thurstone,3,n.obs=213) #we need to specifiy the n.obs for the test to work f2 <- fa(Thurstone,2, n.obs=213) anova(f2,f3) cleanEx() nameEx("bassAckward") ### * bassAckward flush(stderr()); flush(stdout()) ### Name: bassAckward ### Title: The Bass-Ackward factoring algorithm discussed by Goldberg ### Aliases: bassAckward bassAckward.diagram ### Keywords: multivariate models ### ** Examples bassAckward(Thurstone,4,main="Thurstone data set") f.labels <- list(level1=cs(Approach,Avoid),level2=cs(PosAffect,NegAffect,Constraint), level3 = cs(Extraversion,Agreeableness,Neuroticism,Conscientiousness,Openness)) ba <- bassAckward(psychTools::bfi[1:25],c(2,3,5),labels=f.labels, main="bfi data set from psychTools", values=TRUE) print(ba,short=FALSE) #show the items associated with the 5 level solution fa.lookup(ba,dictionary=psychTools::bfi.dictionary) #now show the items associated with the 3 level solution fa.lookup(ba$fa[[2]],dictionary=psychTools::bfi.dictionary) # compare the 3 factor solution to what get by extracting 3 factors directly f3 <- fa(psychTools::bfi[1:25],3) f3$loadings - ba$fa[[2]]$loadings # these are the same #do pca instead of factors just summarize, don't plot summary(bassAckward(psychTools::bfi[1:25],c(1,3,5,7),fm="pca", main="Components",plot=FALSE)) ##not run, but useful example cleanEx() nameEx("best.scales") ### * best.scales flush(stderr()); flush(stdout()) ### Name: bestScales ### Title: A bootstrap aggregation function for choosing most predictive ### unit weighted items ### Aliases: bestItems bestScales BISCUIT biscuit BISCWIT biscwit ### Keywords: models multivariate tree ### ** Examples #This is an example of 'bagging' (bootstrap aggregation) #not run in order to pass the timing tests for Debian at CRAN #bestboot <- bestScales(psychTools::bfi,criteria=cs(gender,age,education), # n.iter=10,dictionary=psychTools::bfi.dictionary[1:3]) #bestboot #compare with 10 fold cross validation cleanEx() nameEx("bfi") ### * bfi flush(stderr()); flush(stdout()) ### Name: bfi ### Title: 25 Personality items representing 5 factors ### Aliases: bfi bfi.keys ### Keywords: datasets ### ** Examples data(bfi) psych::describe(bfi) # create the bfi.keys (actually already saved in the data file) keys <- list(agree=c("-A1","A2","A3","A4","A5"),conscientious=c("C1","C2","C3","-C4","-C5"), extraversion=c("-E1","-E2","E3","E4","E5"),neuroticism=c("N1","N2","N3","N4","N5"), openness = c("O1","-O2","O3","O4","-O5")) scores <- psych::scoreItems(keys,bfi,min=1,max=6) #specify the minimum and maximum values scores #show the use of the fa.lookup with a dictionary #psych::keys.lookup(bfi.keys,bfi.dictionary[,1:4]) #deprecated -- use psychTools cleanEx() nameEx("bi.bars") ### * bi.bars flush(stderr()); flush(stdout()) ### Name: bi.bars ### Title: Draw pairs of bargraphs based on two groups ### Aliases: bi.bars ### Keywords: hplot ### ** Examples #data(bfi) bi.bars(psychTools::bfi,"age","gender" ,ylab="Age",main="Age by males and females") bi.bars(psychTools::bfi,"education","gender",xlab="Education", main="Education by gender",horiz=FALSE) cleanEx() nameEx("bifactor") ### * bifactor flush(stderr()); flush(stdout()) ### Name: Bechtoldt ### Title: Seven data sets showing a bifactor solution. ### Aliases: Bechtoldt.1 Bechtoldt.2 Bechtoldt Holzinger Holzinger.9 Reise ### Thurstone Thurstone.33 Thurstone.9 ### Keywords: datasets ### ** Examples if(!require(GPArotation)) {message("I am sorry, to run omega requires GPArotation") } else { #holz <- omega(Holzinger,4, title = "14 ability tests from Holzinger-Swineford") #bf <- omega(Reise,5,title="16 health items from Reise") #omega(Reise,5,labels=colnames(Reise),title="16 health items from Reise") thur.om <- omega(Thurstone,title="9 variables from Thurstone") #compare with thur.bf <- fa(Thurstone,3,rotate="biquartimin") factor.congruence(thur.om,thur.bf) } cleanEx() nameEx("bigCor") ### * bigCor flush(stderr()); flush(stdout()) ### Name: bigCor ### Title: Find large correlation matrices by stitching together smaller ### ones found more rapidly ### Aliases: bigCor ### Keywords: models multivariate ### ** Examples R <- bigCor(psychTools::bfi,10) #compare the results with r.bfi <- cor(psychTools::bfi,use="pairwise") all.equal(R,r.bfi) cleanEx() nameEx("biplot.psych") ### * biplot.psych flush(stderr()); flush(stdout()) ### Name: biplot.psych ### Title: Draw biplots of factor or component scores by factor or ### component loadings ### Aliases: biplot.psych ### Keywords: multivariate hplot ### ** Examples #the standard example data(USArrests) fa2 <- fa(USArrests,2,scores=TRUE) biplot(fa2,labels=rownames(USArrests)) # plot the 3 factor solution #data(bfi) fa3 <- fa(psychTools::bfi[1:200,1:15],3,scores=TRUE) biplot(fa3) #just plot factors 1 and 3 from that solution biplot(fa3,choose=c(1,3)) # fa2 <- fa(psychTools::bfi[16:25],2) #factor analysis fa2$scores <- fa2$scores[1:100,] #just take the first 100 #now plot with different colors and shapes for males and females biplot(fa2,pch=c(24,21)[psychTools::bfi[1:100,"gender"]], group =psychTools::bfi[1:100,"gender"], main="Biplot of Conscientiousness and Neuroticism by gender") r <- cor(psychTools::bfi[1:200,1:10], use="pairwise") #find the correlations f2 <- fa(r,2) x <- list() x$scores <- factor.scores(psychTools::bfi[1:200,1:10],f2) x$loadings <- f2$loadings class(x) <- c('psych','fa') biplot(x,main="biplot from correlation matrix and factor scores") cleanEx() nameEx("block.random") ### * block.random flush(stderr()); flush(stdout()) ### Name: block.random ### Title: Create a block randomized structure for n independent variables ### Aliases: block.random ### Keywords: multivariate ### ** Examples br <- block.random(n=24,c(2,3)) pairs.panels(br) br <- block.random(96,c(time=4,drug=3,sex=2)) pairs.panels(br) cleanEx() nameEx("bock.table") ### * bock.table flush(stderr()); flush(stdout()) ### Name: bock ### Title: Bock and Liberman (1970) data set of 1000 observations of the ### LSAT ### Aliases: bock bock.table lsat6 lsat7 bock.lsat ### Keywords: datasets ### ** Examples data(bock) responses <- table2df(bock.table[,2:6],count=bock.table[,7], labs= paste("lsat6.",1:5,sep="")) describe(responses) ## maybe str(bock.table) ; plot(bock.table) ... cleanEx() nameEx("cattell") ### * cattell flush(stderr()); flush(stdout()) ### Name: cattell ### Title: 12 cognitive variables from Cattell (1963) ### Aliases: cattell ### Keywords: datasets ### ** Examples data(cattell) corPlot(cattell,numbers=TRUE,upper=FALSE,diag=FALSE, main="12 cognitive variables from Cattell (1963)",xlas=2) 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 <- fa(circ.data,2) plot(circ.fa,title="Circumplex Structure") ct <- circ.tests(circ.fa) #compare with non-circumplex data simp.data <- item.sim(24,500) simp.fa <- fa(simp.data,2) plot(simp.fa,title="Simple Structure") st <- circ.tests(simp.fa) res <- rbind(ct[1:4],st[1:4]) rownames(res) <- c("circumplex","Simple") print(res,digits=2) cleanEx() nameEx("cluster.cor") ### * cluster.cor flush(stderr()); flush(stdout()) ### Name: scoreOverlap ### Title: Find correlations of composite variables (corrected for overlap) ### from a larger matrix. ### Aliases: cluster.cor scoreOverlap scoreBy ### Keywords: multivariate models ### ** Examples #use the msq data set that shows the structure of energetic and tense arousal small.msq <- psychTools::msq[ c("active", "energetic", "vigorous", "wakeful", "wide.awake", "full.of.pep", "lively", "sleepy", "tired", "drowsy","intense", "jittery", "fearful", "tense", "clutched.up", "quiet", "still", "placid", "calm", "at.rest") ] small.R <- cor(small.msq,use="pairwise") keys.list <- list( EA = c("active", "energetic", "vigorous", "wakeful", "wide.awake", "full.of.pep", "lively", "-sleepy", "-tired", "-drowsy"), TA =c("intense", "jittery", "fearful", "tense", "clutched.up", "-quiet", "-still", "-placid", "-calm", "-at.rest") , high.EA = c("active", "energetic", "vigorous", "wakeful", "wide.awake", "full.of.pep", "lively"), low.EA =c("sleepy", "tired", "drowsy"), lowTA= c("quiet", "still", "placid", "calm", "at.rest"), highTA = c("intense", "jittery", "fearful", "tense", "clutched.up") ) keys <- make.keys(small.R,keys.list) adjusted.scales <- scoreOverlap(keys.list,small.R) #compare with unadjusted confounded.scales <- cluster.cor(keys,small.R) summary(adjusted.scales) #note that the EA and high and low EA and TA and high and low TA # scale correlations are confounded summary(confounded.scales) bfi.stats <- statsBy(bfi,group="education",cors=TRUE ,cor="cov") #specify to find covariances bfi.plus.keys <- c(bfi.keys,gender="gender",age ="age") bfi.by <- scoreBy(bfi.plus.keys,bfi.stats) bfi.by$var #to show the variances of each scale by groupl round(bfi.by$cor.mat,2) #the correlations by group bfi.by$alpha 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) fit <- cluster.fit(r.mat,iq.clus$loadings,iq.clus$clusters) fit 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) 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 fa.plot factor.plot ### Keywords: multivariate hplot cluster ### ** Examples circ.data <- circ.sim(24,500) circ.fa <- fa(circ.data,2) plot(circ.fa,cut=.5) f5 <- fa(psychTools::bfi[1:25],5) plot(f5,labels=colnames(psychTools::bfi)[1:25],show.points=FALSE) plot(f5,labels=colnames(psychTools::bfi)[1:25],show.points=FALSE,choose=c(1,2,4)) 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 cluster.cor(keys,test.data) cleanEx() nameEx("cohen.d") ### * cohen.d flush(stderr()); flush(stdout()) ### Name: cohen.d ### Title: Find Cohen d and confidence intervals ### Aliases: cohen.d d.robust cohen.d.ci d.ci cohen.d.by d2r r2d d2t t2d ### m2t m2d d2OVL d2OVL2 d2CL d2U3 cd.validity ### Keywords: models multivariate ### ** Examples cohen.d(sat.act,"gender") #robust version round(d.robust(sat.act,"gender")$robust.d,2) #formula input is nicer cohen.d(sat.act ~ gender) #formula input version #report cohen.d by another group cd <- cohen.d.by(sat.act,"gender","education") cohen.d(SATV + SATQ ~ gender, data=sat.act) #just choose two variables summary(cd) #summarize the output #formula version combines these functions cd <- cohen.d(sat.act ~ gender + education) #find d by gender for each level of education summary(cd) #now show several examples of confidence intervals #one group (d vs 0) #consider the t from the cushny data set t2d( -4.0621,n1=10) d.ci(-1.284549,n1=10) #the confidence interval of the effect of drug on sleep #two groups d.ci(.62,n=64) #equal group size d.ci(.62,n1=35,n2=29) #unequal group size #several examples of d and t from data m2d(52.58,-70.65,49.9,47.5) #Terman and Miles 1936 #graphically show the various overlap statistics curve(d2OVL2(x),0,3,xlab="d",ylab="",lty="dashed", main="Four representations of effect size (d) ") curve(d2OVL(x),0,3,xlab="d",add=TRUE,) curve(d2CL(x),0,3,add=TRUE) curve(d2U3(x), add=TRUE,lty="dotted") text(1,.37,"OVL2") text(2,.37,"OVL") text(1,.88,"U3") text(2, .88,"CL") 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 comorbidity(.2,.15,.1,c("Anxiety","Depression")) cleanEx() nameEx("congruence") ### * congruence flush(stderr()); flush(stdout()) ### Name: congruence ### Title: Matrix and profile congruences and distances ### Aliases: congruence cohen.profile distance ### Keywords: multivariate models ### ** Examples #cohen's example # a and b have reversed one item around the midpoint co <- data.frame(ira=c(2,6,5,6,4), jim=c(1,3,5,4,4), a=c(5,6,5,6,4),b=c(6,3,5,4,4)) lowerMat(congruence(co-3.5)) # congruence around the midpoint is insensitive to reflection lowerCor(co) #the correlation is not lowerMat(congruence(scale(co,scale=FALSE))) #zero centered congruence is r cohen.profile(co) cleanEx() nameEx("cor.ci") ### * cor.ci flush(stderr()); flush(stdout()) ### Name: corCi ### Title: Bootstrapped and normal confidence intervals for raw and ### composite correlations ### Aliases: corCi cor.ci ### Keywords: multivariate models ### ** Examples #find confidence intervals of a correlation matrix with specified sample size ci <- corCi(Thurstone[1:6,1:6],n=213) ci #show them R <- cor.plot.upperLowerCi(ci) #show them graphically R #show them as a matrix #confidence intervals by bootstrapping requires raw data corCi(psychTools::bfi[1:200,1:10]) # just the first 10 variables #The keys have overlapping scales keys <- list(agree=c("-A1","A2","A3","A4","A5"), conscientious= c("C1", "C2","C3","-C4","-C5"),extraversion=c("-E1","-E2","E3","E4","E5"), neuroticism= c("N1", "N2", "N3","N4","N5"), openness = c("O1","-O2","O3","O4","-O5"), alpha=c("-A1","A2","A3","A4","A5","C1","C2","C3","-C4","-C5","N1","N2","N3","N4","N5"), beta = c("-E1","-E2","E3","E4","E5","O1","-O2","O3","O4","-O5") ) #do not correct for item overlap rci <- corCi(psychTools::bfi[1:200,],keys,n.iter=10,main="correlation with overlapping scales") #also shows the graphic -note the overlap #correct for overlap rci <- cor.ci(psychTools::bfi[1:200,],keys,overlap=TRUE, n.iter=10,main="Correct for overlap") #show the confidence intervals ci <- cor.plot.upperLowerCi(rci) #to show the upper and lower confidence intervals ci #print the confidence intervals in matrix form 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 corPlot cor.plot.upperLowerCi corPlotUpperLowerCi ### Keywords: multivariate hplot ### ** Examples corPlot(Thurstone,main="9 cognitive variables from Thurstone") #just blue implies positive manifold #select just some variables to plot corPlot(Thurstone, zlim=c(0,1),main="9 cognitive variables from Thurstone",select=c(1:3,7:9)) #now show a non-symmetric plot corPlot(Thurstone[4:9,1:3], zlim=c(0,1),main="9 cognitive variables from Thurstone",numbers=TRUE,symmetric=FALSE) #Two ways of including stars to show significance #From the raw data corPlot(sat.act,numbers=TRUE,stars=TRUE) #from a correlation matrix with pvals cp <- corr.test(sat.act) #find the correlations and pvals r<- cp$r p <- cp$p corPlot(r,numbers=TRUE,diag=FALSE,stars=TRUE, pval = p,main="Correlation plot with Holm corrected 'significance'") #now red means less than .5 corPlot(mat.sort(Thurstone),TRUE,zlim=c(0,1), main="9 cognitive variables from Thurstone (sorted by factor loading) ") simp <- sim.circ(24) corPlot(cor(simp),main="24 variables in a circumplex") #scale by raw and adjusted probabilities rs <- corr.test(sat.act[1:200,] ) #find the probabilities of the correlations corPlot(r=rs$r,numbers=TRUE,pval=rs$p,main="Correlations scaled by probability values") #Show the upper and lower confidence intervals cor.plot.upperLowerCi(R=rs,numbers=TRUE) #now do this again, but with lighter colors gr <- colorRampPalette(c("#B52127", "white", "#2171B5")) corPlot(r=rs$r,numbers=TRUE,pval=rs$p,main="Correlations scaled by probability values",gr=gr) cor.plot.upperLowerCi(R=rs,numbers=TRUE,gr=gr) #do multiple plots #Also show the xaxis option op <- par(mfrow=c(2,2)) corPlot(psychTools::ability,show.legend=FALSE,keep.par=FALSE,upper=FALSE) f4 <- fa(psychTools::ability,4) corPlot(f4,show.legend=FALSE,keep.par=FALSE,numbers=TRUE,xlas=3) om <- omega(psychTools::ability,4) corPlot(om,show.legend=FALSE,keep.par=FALSE,numbers=TRUE,xaxis=3) par(op) corPlotUpperLowerCi(rs,adjust=TRUE,main="Holm adjusted confidence intervals",gr=gr) graphics::par(get("par.postscript", pos = 'CheckExEnv')) cleanEx() nameEx("cor.smooth") ### * cor.smooth flush(stderr()); flush(stdout()) ### Name: cor.smooth ### Title: Smooth a non-positive definite correlation matrix to make it ### positive definite ### Aliases: cor.smooth cor.smoother ### Keywords: multivariate models ### ** Examples burt <- psychTools::burt bs <- cor.smooth(psychTools::burt) #burt data set is not positive definite plot(burt[lower.tri(burt)],bs[lower.tri(bs)],ylab="smoothed values",xlab="original values") abline(0,1,lty="dashed") round(burt - bs,3) fa(burt,2) #this throws a warning that the matrix yields an improper solution #Smoothing first throws a warning that the matrix was improper, #but produces a better solution fa(cor.smooth(burt),2) #this next example is a correlation matrix from DeLeuw used as an example #in Knol and ten Berge. #the example is also used in the nearcor documentation cat("pr is the example matrix used in Knol DL, ten Berge (1989)\n") pr <- matrix(c(1, 0.477, 0.644, 0.478, 0.651, 0.826, 0.477, 1, 0.516, 0.233, 0.682, 0.75, 0.644, 0.516, 1, 0.599, 0.581, 0.742, 0.478, 0.233, 0.599, 1, 0.741, 0.8, 0.651, 0.682, 0.581, 0.741, 1, 0.798, 0.826, 0.75, 0.742, 0.8, 0.798, 1), nrow = 6, ncol = 6) sm <- cor.smooth(pr) resid <- pr - sm # several goodness of fit tests # from Knol and ten Berge tr(resid %*% t(resid)) /2 # from nearPD sum(resid^2)/2 cleanEx() nameEx("cor.wt") ### * cor.wt flush(stderr()); flush(stdout()) ### Name: cor.wt ### Title: The sample size weighted correlation may be used in correlating ### aggregated data ### Aliases: cor.wt ### Keywords: multivariate ### ** Examples means.by.age <- statsBy(sat.act,"age") wt.cors <- cor.wt(means.by.age) lowerMat(wt.cors$r) #show the weighted correlations unwt <- lowerCor(means.by.age$mean) mixed <- lowerUpper(unwt,wt.cors$r) #combine both results cor.plot(mixed,TRUE,main="weighted versus unweighted correlations") diff <- lowerUpper(unwt,wt.cors$r,TRUE) cor.plot(diff,TRUE,main="differences of weighted versus unweighted correlations") cleanEx() nameEx("corFiml") ### * corFiml flush(stderr()); flush(stdout()) ### Name: corFiml ### Title: Find a Full Information Maximum Likelihood (FIML) correlation or ### covariance matrix from a data matrix with missing data ### Aliases: corFiml ### Keywords: multivariate models ### ** Examples rML <- corFiml(psychTools::bfi[20:27]) rpw <- cor(psychTools::bfi[20:27],use="pairwise") round(rML - rpw,3) mp <- corFiml(psychTools::bfi[20:27],show=TRUE) mp 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 corr.p ### Keywords: multivariate models ### ** Examples ct <- corr.test(attitude) #find the correlations and give the probabilities ct #show the results cts <- corr.test(attitude[1:3],attitude[4:6]) #reports all values corrected for multiple tests #corr.test(sat.act[1:3],sat.act[4:6],adjust="none") #don't adjust the probabilities #take correlations and show the probabilities as well as the confidence intervals print(corr.p(cts$r,n=30),short=FALSE) #don't adjust the probabilities print(corr.test(sat.act[1:3],sat.act[4:6],adjust="none"),short=FALSE) #print out the stars object without showing quotes print(corr.test(attitude)$stars,quote=FALSE) #note that the adjusted ps are given as well kendall.r <- corr.test(bfi[1:40,4:6], method="kendall", normal=FALSE) #compare with cor.test(x=bfi[1:40,4],y=bfi[1:40,6],method="kendall", exact=FALSE) print(kendall.r,digits=6) 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) # #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)) # or clusters <- matrix(c(rep(1,3),rep(0,7),rep(1,4)),ncol=2) cluster.cor(clusters,cor(attitude)) # #best keys <- make.keys(attitude,list(first=1:3,second=4:7)) scores <- scoreItems(keys,attitude) scores$corrected #However, to do the more general case of correcting correlations for reliabilty #corrected <- cor2cov(x.cor,1/alpha) #diag(corrected) <- 1 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 #data(bfi) cortest.bartlett(psychTools::bfi[1:200,1:10]) #not an identity matrix f3 <- fa(Thurstone,3) f3r <- f3$resid cortest.bartlett(f3r,n=213,diag=FALSE) #incorrect cortest.bartlett(f3r,n=213,diag=TRUE) #correct (by default) 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) cortest.normal(x) #just test if this matrix is an identity x <- sim.congeneric(loads =c(.9,.8,.7,.6,.5),N=1000,short=FALSE) y <- sim.congeneric(loads =c(.9,.8,.7,.6,.5),N=1000,short=FALSE) cortest.normal(x$r,y$r,n1=1000,n2=1000) #The Steiger test cortest.jennrich(x$r,y$r,n1=100,n2=1000) # The Jennrich test cortest.mat(x$r,y$r,n1=1000,n2=1000) #twice the degrees of freedom as the Jennrich cleanEx() nameEx("cosinor") ### * cosinor flush(stderr()); flush(stdout()) ### Name: cosinor ### Title: Functions for analysis of circadian or diurnal data ### Aliases: cosinor circadian.phase cosinor.plot cosinor.period ### circadian.mean circadian.sd circadian.cor circadian.linear.cor ### circadian.stats circadian.F circadian.reliability circular.mean ### circular.cor ### Keywords: multivariate ### ** Examples time <- seq(1:24) #create a 24 hour time pure <- matrix(time,24,18) colnames(pure) <- paste0("H",1:18) pure <- data.frame(time,cos((pure - col(pure))*pi/12)*3 + 3) #18 different phases but scaled to 0-6 match mood data matplot(pure[-1],type="l",main="Pure circadian arousal rhythms", xlab="time of day",ylab="Arousal") op <- par(mfrow=c(2,2)) cosinor.plot(1,3,pure) cosinor.plot(1,5,pure) cosinor.plot(1,8,pure) cosinor.plot(1,12,pure) p <- cosinor(pure) #find the acrophases (should match the input) #now, test finding the acrophases for different subjects on 3 variables #They should be the first 3, second 3, etc. acrophases of pure pp <- matrix(NA,nrow=6*24,ncol=4) pure <- as.matrix(pure) pp[,1] <- rep(pure[,1],6) pp[1:24,2:4] <- pure[1:24,2:4] pp[25:48,2:4] <- pure[1:24,5:7] *2 #to test different variances pp[49:72,2:4] <- pure[1:24,8:10] *3 pp[73:96,2:4] <- pure[1:24,11:13] pp[97:120,2:4] <- pure[1:24,14:16] pp[121:144,2:4] <- pure[1:24,17:19] pure.df <- data.frame(ID = rep(1:6,each=24),pp) colnames(pure.df) <- c("ID","Time",paste0("V",1:3)) cosinor("Time",3:5,"ID",pure.df) op <- par(mfrow=c(2,2)) cosinor.plot(2,3,pure.df,IDloc=1,ID="1") cosinor.plot(2,3,pure.df,IDloc=1,ID="2") cosinor.plot(2,3,pure.df,IDloc=1,ID="3") cosinor.plot(2,3,pure.df,IDloc=1,ID="4") #now, show those in one panel as spagetti plots op <- par(mfrow=c(1,1)) cosinor.plot(2,3,pure.df,IDloc=1,ID="1",multi=TRUE,ylim=c(0,20),ylab="Modeled") cosinor.plot(2,3,pure.df,IDloc=1,ID="2",multi=TRUE,add=TRUE,lty="dotdash") cosinor.plot(2,3,pure.df,IDloc=1,ID="3",multi=TRUE,add=TRUE,lty="dashed") cosinor.plot(2,3,pure.df,IDloc=1,ID="4",multi=TRUE,add=TRUE,lty="dotted") set.seed(42) #what else? noisy <- pure noisy[,2:19]<- noisy[,2:19] + rnorm(24*18,0,.2) n <- cosinor(time,noisy) #add a bit of noise small.pure <- pure[c(8,11,14,17,20,23),] small.noisy <- noisy[c(8,11,14,17,20,23),] small.time <- c(8,11,14,17,20,23) cosinor.plot(1,3,small.pure,multi=TRUE) cosinor.plot(1,3,small.noisy,multi=TRUE,add=TRUE,lty="dashed") # sp <- cosinor(small.pure) # spo <- cosinor(small.pure,opti=TRUE) #iterative fit # sn <- cosinor(small.noisy) #linear # sno <- cosinor(small.noisy,opti=TRUE) #iterative # sum.df <- data.frame(pure=p,noisy = n, small=sp,small.noise = sn, # small.opt=spo,small.noise.opt=sno) # round(sum.df,2) # round(circadian.cor(sum.df[,c(1,3,5,7,9,11)]),2) #compare alternatives # # #now, lets form three "subjects" and show how the grouping variable works # mixed.df <- rbind(small.pure,small.noisy,noisy) # mixed.df <- data.frame(ID=c(rep(1,6),rep(2,6),rep(3,24)), # time=c(rep(c(8,11,14,17,20,23),2),1:24),mixed.df) # group.df <- cosinor(angle="time",x=2:20,code="ID",data=mixed.df) # round(group.df,2) #compare these values to the sp,sn,and n values done separately graphics::par(get("par.postscript", pos = 'CheckExEnv')) cleanEx() nameEx("count.pairwise") ### * count.pairwise flush(stderr()); flush(stdout()) ### Name: pairwiseCount ### Title: Count number of pairwise cases for a data set with missing (NA) ### data and impute values. ### Aliases: pairwiseCount pairwiseCountBig count.pairwise pairwiseDescribe ### pairwiseZero pairwiseSample pairwiseReport pairwiseImpute ### pairwisePlot ### Keywords: models multivariate ### ** Examples x <- matrix(rnorm(900),ncol=6) y <- matrix(rnorm(450),ncol=3) x[x < 0] <- NA y[y > 1] <- NA pairwiseCount(x) pairwiseCount(y) pairwiseCount(x,y) pairwiseCount(x,diagonal=FALSE) pairwiseDescribe(x,quant=c(.1,.25,.5,.75,.9)) #examine the structure of the ability data set keys <- list(ICAR16=colnames(psychTools::ability),reasoning = cs(reason.4,reason.16,reason.17,reason.19), letters=cs(letter.7, letter.33,letter.34,letter.58, letter.7), matrix=cs(matrix.45,matrix.46,matrix.47,matrix.55), rotate=cs(rotate.3,rotate.4,rotate.6,rotate.8)) pairwiseImpute(keys,psychTools::ability) cleanEx() nameEx("cta") ### * cta flush(stderr()); flush(stdout()) ### Name: cta ### Title: Simulate the C(ues) T(endency) A(ction) model of motivation ### Aliases: cta cta.15 ### 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 #these next are examples without graphic output #not run #two introverts #c2i <- c(.95,1.05) #cta(n=2,t=10000,cues=c2i,type="none") #two extraverts #c2e <- c(3.95,4.05) #cta(n=2,t=10000,cues=c2e,type="none") #three introverts #c3i <- c(.95,1,1.05) #cta(3,t=10000,cues=c3i,type="none") #three extraverts #c3i <- c(3.95,4, 4.05) #cta(3,10000,c3i,type="none") #mixed #c3 <- c(1,2.5,4) #cta(3,10000,c3,type="none") cleanEx() nameEx("densityBy") ### * densityBy flush(stderr()); flush(stdout()) ### Name: densityBy ### Title: Create a 'violin plot' or density plot of the distribution of a ### set of variables ### Aliases: densityBy violinBy violin ### Keywords: multivariate hplot ### ** Examples violin(psychTools::bfi[4:8]) violin(SATV + SATQ ~ gender, data=sat.act, grp.name =cs(MV,FV,MQ,FQ)) #formula input violinBy(psychTools::bfi,var=4:7,grp ="gender",grp.name=c("M","F")) densityBy(SATV ~ gender,data =sat.act,legend=1) #formula input cleanEx() nameEx("deprecated") ### * deprecated flush(stderr()); flush(stdout()) ### Name: fa.poly ### Title: Deprecated Exploratory Factor analysis functions. Please use fa ### Aliases: factor.pa factor.minres factor.wls fa.poly ### Keywords: multivariate models ### ** Examples #none, you should see fa #using the Harman 24 mental tests, compare a principal factor with a principal components solution cleanEx() nameEx("describe") ### * describe flush(stderr()); flush(stdout()) ### Name: describe ### Title: Basic descriptive statistics useful for psychometrics ### Aliases: describe describeData describeFast ### Keywords: multivariate models univar ### ** Examples data(sat.act) describe(sat.act) describe(sat.act ~ gender) #formula mode option calls describeBy for the entire data frame describe(SATV + SATQ ~ gender, data=sat.act) #formula mode specifies just two variables describe(sat.act,skew=FALSE) describe(sat.act,IQR=TRUE) #show the interquartile Range describe(sat.act,quant=c(.1,.25,.5,.75,.90) ) #find the 10th, 25th, 50th, #75th and 90th percentiles describeData(sat.act) #the fast version just gives counts and head and tail print(describeFast(sat.act),short=FALSE) #even faster is just counts (just less information) #now show how to adjust the displayed number of digits des <- describe(sat.act) #find the descriptive statistics. Keep the original accuracy des #show the normal output, which is rounded to 2 decimals print(des,digits=3) #show the output, but round to 3 (trailing) digits print(des, signif=3) #round all numbers to the 3 significant digits cleanEx() nameEx("describe.by") ### * describe.by flush(stderr()); flush(stdout()) ### Name: describeBy ### Title: Basic summary statistics by group ### Aliases: describeBy describe.by ### Keywords: models univar ### ** Examples data(sat.act) describeBy(sat.act,sat.act$gender) #just one grouping variable describeBy(sat.act ~ gender) #describe the entire set formula input describeBy(SATV + SATQ ~ gender,data =sat.act) #specify the data set if using formula #describeBy(sat.act,list(sat.act$gender,sat.act$education)) #two grouping variables describeBy(sat.act ~ gender + education) #two grouping variables des.mat <- describeBy(age ~ education,mat=TRUE,data = sat.act) #matrix (data.frame) output des.mat <- describeBy(age ~ education + gender, data=sat.act, mat=TRUE,digits=2) #matrix output rounded to 2 decimals cleanEx() nameEx("diagram") ### * diagram flush(stderr()); flush(stdout()) ### Name: diagram ### Title: Helper functions for drawing path model diagrams ### Aliases: diagram dia.rect dia.ellipse dia.ellipse1 dia.arrow dia.curve ### dia.curved.arrow dia.self dia.shape dia.triangle dia.cone multi.rect ### multi.arrow multi.curved.arrow multi.self ### Keywords: multivariate hplot ### ** Examples #first, show the primitives xlim=c(-2,10) ylim=c(0,10) plot(NA,xlim=xlim,ylim=ylim,main="Demonstration of diagram functions",axes=FALSE,xlab="",ylab="") ul <- dia.rect(1,9,labels="upper left",xlim=xlim,ylim=ylim) ml <- dia.rect(1,6,"middle left",xlim=xlim,ylim=ylim) ll <- dia.rect(1,3,labels="lower left",xlim=xlim,ylim=ylim) bl <- dia.rect(1,1,"bottom left",xlim=xlim,ylim=ylim) lr <- dia.ellipse(7,3,"lower right",xlim=xlim,ylim=ylim,e.size=.07) ur <- dia.ellipse(7,9,"upper right",xlim=xlim,ylim=ylim,e.size=.07) mr <- dia.ellipse(7,6,"middle right",xlim=xlim,ylim=ylim,e.size=.07) lm <- dia.triangle(4,1,"Lower Middle",xlim=xlim,ylim=ylim) br <- dia.rect(9,1,"bottom right",xlim=xlim,ylim=ylim) dia.curve(from=ul$left,to=bl$left,"double headed",scale=-1) dia.arrow(from=lr,to=ul,labels="right to left") dia.arrow(from=ul,to=ur,labels="left to right") dia.curved.arrow(from=lr,to=ll,labels ="right to left") dia.curved.arrow(to=ur,from=ul,labels ="left to right") dia.curve(ll$top,ul$bottom,"right") #for rectangles, specify where to point dia.curve(ll$top,ul$bottom,"left",scale=-1) #for rectangles, specify where to point dia.curve(mr,ur,"up") #but for ellipses, you may just point to it. dia.curve(mr,lr,"down") dia.curve(mr,ur,"up") dia.curved.arrow(mr,ur,"up") #but for ellipses, you may just point to it. dia.curved.arrow(mr,lr,"down") #but for ellipses, you may just point to it. dia.curved.arrow(ur$right,mr$right,"3") dia.curve(ml,mr,"across") dia.curve(ur$right,lr$right,"top down",scale =2) dia.curved.arrow(br$top,lr$right,"up") dia.curved.arrow(bl,br,"left to right") dia.curved.arrow(br,bl,"right to left",scale=-1) dia.arrow(bl,ll$bottom) dia.curved.arrow(ml,ll$right) dia.curved.arrow(mr,lr$top) #now, put them together in a factor analysis diagram v9 <- sim.hierarchical() f3 <- fa(v9,3,rotate="cluster") fa.diagram(f3,error=TRUE,side=3) cleanEx() nameEx("draw.tetra") ### * draw.tetra flush(stderr()); flush(stdout()) ### Name: draw.tetra ### Title: Draw a correlation ellipse and two normal curves to demonstrate ### tetrachoric correlation ### Aliases: draw.tetra draw.cor ### Keywords: multivariate hplot ### ** Examples #if(require(mvtnorm)) { #draw.tetra(.5,1,1) #draw.tetra(.8,2,1)} else {print("draw.tetra requires the mvtnorm package") #draw.cor(.5,cuts=c(0,0))} draw.tetra(.5,1,1) draw.tetra(.8,2,1) draw.cor(.5,cuts=c(0,0)) cleanEx() nameEx("dummy.code") ### * dummy.code flush(stderr()); flush(stdout()) ### Name: dummy.code ### Title: Create dummy coded variables ### Aliases: dummy.code ### Keywords: multivariate models ### ** Examples new <- dummy.code(sat.act$education) new.sat <- data.frame(new,sat.act) round(cor(new.sat,use="pairwise"),2) #dum.smoke <- dummy.code(spi$smoke,group=2:9) #table(dum.smoke,spi$smoke) #dum.age <- dummy.code(round(spi$age/5)*5,top=5) #the most frequent five year blocks cleanEx() nameEx("dwyer") ### * dwyer flush(stderr()); flush(stdout()) ### Name: Dwyer ### Title: 8 cognitive variables used by Dwyer for an example. ### Aliases: Dwyer ### Keywords: datasets ### ** Examples data(Dwyer) Ro <- Dwyer[1:7,1:7] Roe <- Dwyer[1:7,8] fo <- fa(Ro,2,rotate="none") fa.extension(Roe,fo) 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 y <- princomp(covmat=Harman74.cor$cov) y$loadings[1:8,1:4] #as they appear from princomp eigen.loadings(x)[1:8,1:4] # rescaled by the eigen values cleanEx() nameEx("ellipses") ### * ellipses flush(stderr()); flush(stdout()) ### Name: ellipses ### Title: Plot data and 1 and 2 sigma correlation ellipses ### Aliases: ellipses minkowski ### Keywords: multivariate hplot ### ** Examples data(psychTools::galton) galton <- psychTools::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 minkowski(2,main="Minkowski circles") minkowski(1,TRUE) minkowski(4,TRUE) cleanEx() nameEx("error.bars") ### * error.bars flush(stderr()); flush(stdout()) ### Name: error.bars ### Title: Plot means and confidence intervals ### Aliases: error.bars error.bars.tab ### Keywords: multivariate hplot ### ** Examples set.seed(42) x <- matrix(rnorm(1000),ncol=20) boxplot(x,notch=TRUE,main="Notched boxplot with error bars") error.bars(x,add=TRUE) abline(h=0) #show 50% confidence regions and color each variable separately error.bars(attitude,alpha=.5, main="50 percent confidence limits",col=rainbow(ncol(attitude)) ) error.bars(attitude,bar=TRUE) #show the use of bar graphs #combine with a strip chart and boxplot stripchart(attitude,vertical=TRUE,method="jitter",jitter=.1,pch=19, main="Stripchart with 95 percent confidence limits") boxplot(attitude,add=TRUE) error.bars(attitude,add=TRUE,arrow.len=.2) #use statistics from somewhere else #by specifying n, we are using the t distribution for confidences #The first example allows the variables to be spaced along the x axis my.stats <- data.frame(values=c(1,2,8),mean=c(10,12,18),se=c(2,3,5),n=c(5,10,20)) error.bars(stats=my.stats,type="b",main="data with confidence intervals") #don't connect the groups my.stats <- data.frame(values=c(1,2,8),mean=c(10,12,18),se=c(2,3,5),n=c(5,10,20)) error.bars(stats=my.stats,main="data with confidence intervals") #by not specifying value, the groups are equally spaced my.stats <- data.frame(mean=c(10,12,18),se=c(2,3,5),n=c(5,10,20)) rownames(my.stats) <- c("First", "Second","Third") error.bars(stats=my.stats,xlab="Condition",ylab="Score") #Consider the case where we get stats from describe temp <- describe(attitude) error.bars(stats=temp) #show these do not differ from the other way by overlaying the two error.bars(attitude,add=TRUE,col="red") #n is omitted #the error distribution is a normal distribution my.stats <- data.frame(mean=c(2,4,8),se=c(2,1,2)) rownames(my.stats) <- c("First", "Second","Third") error.bars(stats=my.stats,xlab="Condition",ylab="Score") #n is specified #compare this with small n which shows larger confidence regions my.stats <- data.frame(mean=c(2,4,8),se=c(2,1,2),n=c(10,10,3)) rownames(my.stats) <- c("First", "Second","Third") error.bars(stats=my.stats,xlab="Condition",ylab="Score") #example of arrest rates (as percentage of condition) arrest <- data.frame(Control=c(14,21),Treated =c(3,23)) rownames(arrest) <- c("Arrested","Not Arrested") error.bars.tab(arrest,ylab="Probability of Arrest",xlab="Control vs Treatment", main="Probability of Arrest varies by treatment") #Show the raw rates error.bars.tab(arrest,raw=TRUE,ylab="Number Arrested",xlab="Control vs Treatment", main="Count of Arrest varies by treatment") #If a grouping variable is specified, the function calls error.bars.by #Use error.bars.by to have more control over the output. #Show how to use grouping variables error.bars(SATV + SATQ ~ gender, data=sat.act) #one grouping variable, formula input error.bars(SATV + SATQ ~ education + gender,data=sat.act)#two grouping variables 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) #The generic plot of variables by group error.bars.by( SATV + SATQ ~ gender,data=sat.act) #formula input error.bars.by( SATV + SATQ ~ gender,data=sat.act,v.lab=cs(male,female)) #labels error.bars.by(SATV + SATQ ~ education + gender, data =sat.act) #see below error.bars.by(sat.act[1:4],sat.act$gender,legend=7) #specification of variables error.bars.by(sat.act[1:4],sat.act$gender,legend=7,labels=cs(male,female)) #a bar plot error.bars.by(sat.act[5:6],sat.act$gender,bars=TRUE,labels=c("male","female"), main="SAT V and SAT Q by gender",ylim=c(0,800),colors=c("red","blue"), legend=5,v.labels=c("SATV","SATQ")) #draw a barplot #a bar plot of SAT by age -- not recommended, see the next plot error.bars.by(SATV + SATQ ~ education,data=sat.act,bars=TRUE,xlab="Education", main="95 percent confidence limits of Sat V and Sat Q", ylim=c(0,800), v.labels=c("SATV","SATQ"),colors=c("red","blue") ) #a better graph uses points not bars #use formulat input #plot SAT V and SAT Q by education error.bars.by(SATV + SATQ ~ education,data=sat.act,TRUE, xlab="Education", legend=5,labels=colnames(sat.act[5:6]),ylim=c(525,700), main="self reported SAT scores by education", v.lab =c("HS","in coll", "< 16", "BA/BS", "in Grad", "Grad/Prof")) #make the cats eyes semi-transparent by specifying a negative density error.bars.by(SATV + SATQ ~ education,data=sat.act, xlab="Education", legend=5,labels=c("SATV","SATQ"),ylim=c(525,700), main="self reported SAT scores by education",density=-10, v.lab =c("HS","in coll", "< 16", "BA/BS", "in Grad", "Grad/Prof")) #use labels to specify the 2nd grouping variable, v.lab to specify the first error.bars.by(SATV ~ education + gender,data=sat.act, xlab="Education", legend=5,labels=cs(male,female),ylim=c(525,700), main="self reported SAT scores by education",density=-10, v.lab =c("HS","in coll", "< 16", "BA/BS", "in Grad", "Grad/Prof"), colors=c("red","blue")) #now for a more complicated examples using 25 big 5 items scored into 5 scales #and showing age trends by decade #this shows how to convert many levels of a grouping variable (age) into more manageable levels. data(bfi) #The Big 5 data #first create the keys 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(psychTools::bfi,keys.list) #then create the scores for those older than 10 and less than 80 bfis <- subset(psychTools::bfi,((psychTools::bfi$age > 10) & (psychTools::bfi$age < 80))) scores <- scoreItems(keys,bfis,min=1,max=6) #set the right limits for item reversals #now draw the results by age #specify the particular colors to use error.bars.by(scores$scores,round(bfis$age/10)*10,by.var=TRUE, main="BFI age trends",legend=3,labels=colnames(scores$scores), xlab="Age",ylab="Mean item score", colors=cs(green,yellow,black,red,blue), v.labels =cs(10-14,15-24,25-34,35-44,45-54,55-64,65-74)) #show transparency error.bars.by(scores$scores,round(bfis$age/10)*10,by.var=TRUE, main="BFI age trends",legend=3,labels=colnames(scores$scores), xlab="Age",ylab="Mean item score", density=-10, colors=cs(green,yellow,black,red,blue), v.labels =cs(10-14,15-24,25-34,35-44,45-54,55-64,65-74)) cleanEx() nameEx("error.circles") ### * error.circles flush(stderr()); flush(stdout()) ### Name: errorCircles ### Title: Two way plots of means, error bars, and sample sizes ### Aliases: errorCircles ### Keywords: multivariate hplot ### ** Examples #BFI scores for males and females errorCircles(1:25,1:25,data=psychTools::bfi,group="gender",paired=TRUE,ylab="female scores", xlab="male scores",main="BFI scores by gender") abline(a=0,b=1) #drop the circles since all samples are the same sizes errorCircles(1:25,1:25,data=psychTools::bfi,group="gender",paired=TRUE,circles=FALSE, ylab="female scores",xlab="male scores",main="BFI scores by gender") abline(a=0,b=1) data(psychTools::affect) colors <- c("black","red","white","blue") films <- c("Sad","Horror","Neutral","Happy") affect.stats <- errorCircles("EA2","TA2",data=psychTools::affect[-c(1,20)], group="Film",labels=films, xlab="Energetic Arousal",ylab="Tense Arousal",ylim=c(10,22),xlim=c(8,20), pch=16,cex=2,colors=colors, main ="EA and TA pre and post affective movies") #now, use the stats from the prior run errorCircles("EA1","TA1",data=affect.stats,labels=films,pch=16,cex=2,colors=colors,add=TRUE) #show sample size with the size of the circles errorCircles("SATV","SATQ",sat.act,group="education") #Can also provide error.bars.by functionality errorCircles(2,5,group=2,data=sat.act,circles=FALSE,pch=16,colors="blue", ylim= c(200,800),main="SATV by education",labels="") #just do the breakdown and then show the points # errorCircles(3,5,group=3,data=sat.act,circles=FALSE,pch=16,colors="blue", # ylim= c(200,800),main="SATV by age",labels="",bars=FALSE) 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 #just draw one pair of variables desc <- describe(attitude) x <- desc[1,] y <- desc[2,] error.crosses(x,y,xlab=rownames(x),ylab=rownames(y)) #now for a bit more complicated plotting data(psychTools::bfi) desc <- describeBy(psychTools::bfi[1:25],psychTools::bfi$gender) #select a high and low group error.crosses(desc$'1',desc$'2',ylab="female scores", xlab="male scores",main="BFI scores by gender") abline(a=0,b=1) #do it from summary statistics (using standard errors) g1.stats <- data.frame(n=c(10,20,30),mean=c(10,12,18),se=c(2,3,5)) g2.stats <- data.frame(n=c(15,20,25),mean=c(6,14,15),se =c(1,2,3)) error.crosses(g1.stats,g2.stats) #Or, if you prefer to draw +/- 1 sd. instead of 95% confidence g1.stats <- data.frame(n=c(10,20,30),mean=c(10,12,18),sd=c(2,3,5)) g2.stats <- data.frame(n=c(15,20,25),mean=c(6,14,15),sd =c(1,2,3)) error.crosses(g1.stats,g2.stats,sd=TRUE) #and seem even fancy plotting: This is taken from a study of mood #four films were given (sad, horror, neutral, happy) #with a pre and post test data(psychTools::affect) colors <- c("black","red","green","blue") films <- c("Sad","Horror","Neutral","Happy") affect.mat <- describeBy(psychTools::affect[10:17],psychTools::affect$Film,mat=TRUE) error.crosses(affect.mat[c(1:4,17:20),],affect.mat[c(5:8,21:24),], labels=films[affect.mat$group1],xlab="Energetic Arousal", ylab="Tense Arousal",colors = colors[affect.mat$group1],pch=16,cex=2) cleanEx() nameEx("error.dots") ### * error.dots flush(stderr()); flush(stdout()) ### Name: error.dots ### Title: Show a dot.chart with error bars for different groups or ### variables ### Aliases: error.dots ### Keywords: multivariate hplot ### ** Examples temp <- error.dots(psychTools::bfi[1:25],sort=TRUE, xlab="Mean score for the item, sorted by difficulty") error.dots(psychTools::bfi[1:25],sort=TRUE, order=temp$order, add=TRUE, eyes=TRUE) #over plot with eyes error.dots(psychTools::ability,eyes=TRUE, xlab="Mean score for the item") cd <- cohen.d(psychTools::bfi[1:26],"gender") temp <- error.dots(cd, select=c(1:15,21:25),head=12,tail=13, main="Cohen d and confidence intervals of BFI by gender") error.dots(cd,select=c(16:20),head=13,tail=12,col="blue",add=TRUE,fg="red" ,main="") abline(v=0) cleanEx() nameEx("esem") ### * esem flush(stderr()); flush(stdout()) ### Name: esem ### Title: Perform and Exploratory Structural Equation Model (ESEM) by ### using factor extension techniques ### Aliases: esem esem.diagram interbattery ### Keywords: multivariate models ### ** Examples #make up a sem like problem using sim.structure fx <-matrix(c( .9,.8,.6,rep(0,4),.6,.8,-.7),ncol=2) fy <- matrix(c(.6,.5,.4),ncol=1) rownames(fx) <- c("V","Q","A","nach","Anx") rownames(fy)<- c("gpa","Pre","MA") Phi <-matrix( c(1,0,.7,.0,1,.7,.7,.7,1),ncol=3) gre.gpa <- sim.structural(fx,Phi,fy) print(gre.gpa) #now esem it: example <- esem(gre.gpa$model,varsX=1:5,varsY=6:8,nfX=2,nfY=1,n.obs=1000,plot=FALSE) example esem.diagram(example,simple=FALSE) #compare two alternative solutions to the first 2 factors of the neo. #solution 1 is the normal 2 factor solution. #solution 2 is an esem with 1 factor for the first 6 variables, and 1 for the second 6. f2 <- fa(psychTools::neo[1:12,1:12],2) es2 <- esem(psychTools::neo,1:6,7:12,1,1) summary(f2) summary(es2) fa.congruence(f2,es2) interbattery(Thurstone.9,1:4,5:9,2,2) #compare to the solution of Tucker. We are not there yet. cleanEx() nameEx("fa") ### * fa flush(stderr()); flush(stdout()) ### Name: fa ### Title: Exploratory Factor analysis using MinRes (minimum residual) as ### well as EFA by Principal Axis, Weighted Least Squares or Maximum ### Likelihood ### Aliases: fa fac fa.sapa fa.pooled ### 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") #principal components pa <- fa(Harman74.cor$cov,4,fm="pa" ,rotate="varimax") #principal axis uls <- fa(Harman74.cor$cov,4,rotate="varimax") #unweighted least squares is minres wls <- fa(Harman74.cor$cov,4,fm="wls") #weighted least squares #to show the loadings sorted by absolute value print(uls,sort=TRUE) #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)) #note that the order of factors and the sign of some of factors may differ #finally, compare the unrotated factor, ml, uls, and wls solutions wls <- fa(Harman74.cor$cov,4,rotate="none",fm="wls") pa <- fa(Harman74.cor$cov,4,rotate="none",fm="pa") minres <- factanal(factors=4,covmat=Harman74.cor$cov,rotation="none") mle <- fa(Harman74.cor$cov,4,rotate="none",fm="mle") uls <- fa(Harman74.cor$cov,4,rotate="none",fm="uls") factor.congruence(list(minres,mle,pa,wls,uls)) #in particular, note the similarity of the mle and min res solutions #note that the order of factors and the sign of some of factors may 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 Thurstone.33 <- as.matrix(Thurstone.33) mle2 <- fa(Thurstone.33,2,rotate="none",fm="mle") mle3 <- fa(Thurstone.33,3 ,rotate="none",fm="mle") pa2 <- fa(Thurstone.33,2,rotate="none",fm="pa") pa3 <- fa(Thurstone.33,3,rotate="none",fm="pa") mr2 <- fa(Thurstone.33,2,rotate="none") mr3 <- fa(Thurstone.33,3,rotate="none") factor.congruence(list(mle2,mr2,pa2,mle3,pa3,mr3)) #f5 <- fa(psychTools::bfi[1:25],5) #f5 #names are not in ascending numerical order (see note) #colnames(f5$loadings) <- paste("F",1:5,sep="") #f5 #Get the variance accounted for object from the print function p <- print(mr3) round(p$Vaccounted,2) #pool data and fa the pooled result (not run) #datasets.list <- list(bfi[1:500,1:25],bfi[501:1000,1:25], # bfi[1001:1500,1:25],bfi[1501:2000,1:25],bfi[2001:2500,1:25]) #five different data sets #temp <- fa.pooled(datasets.list,5) #do 5 factor analyses, pool the results cleanEx() nameEx("fa.diagram") ### * fa.diagram flush(stderr()); flush(stdout()) ### Name: fa.diagram ### Title: Graph factor loading matrices ### Aliases: fa.graph fa.rgraph fa.diagram extension.diagram het.diagram ### Keywords: multivariate hplot ### ** Examples test.simple <- fa(item.sim(16),2,rotate="oblimin") #if(require(Rgraphviz)) {fa.graph(test.simple) } fa.diagram(test.simple) f3 <- fa(Thurstone,3,rotate="cluster") fa.diagram(f3,cut=.4,digits=2) f3l <- f3$loadings fa.diagram(f3l,main="input from a matrix") Phi <- f3$Phi fa.diagram(f3l,Phi=Phi,main="Input from a matrix") fa.diagram(ICLUST(Thurstone,2,title="Two cluster solution of Thurstone"),main="Input from ICLUST") het.diagram(Thurstone,levels=list(1:4,5:8,3:7)) cleanEx() nameEx("fa.extension") ### * fa.extension flush(stderr()); flush(stdout()) ### Name: fa.extension ### Title: Apply Dwyer's factor extension to find factor loadings for ### extended variables ### Aliases: fa.extension fa.extend ### Keywords: multivariate ### ** Examples #The Dwyer Example Ro <- Dwyer[1:7,1:7] Roe <- Dwyer[1:7,8] fo <- fa(Ro,2,rotate="none") fe <- fa.extension(Roe,fo) #an example from simulated data set.seed(42) d <- sim.item(12) #two orthogonal factors R <- cor(d) Ro <- R[c(1,2,4,5,7,8,10,11),c(1,2,4,5,7,8,10,11)] Roe <- R[c(1,2,4,5,7,8,10,11),c(3,6,9,12)] fo <- fa(Ro,2) fe <- fa.extension(Roe,fo) fa.diagram(fo,fe=fe) #alternatively just specify the original variables and the extension variables fe = fa.extend(R, 2, ov =c(1,2,4,5,7,8,10,11), ev=c(3,6,9,12)) fa.diagram(fe$fo, fe = fe$fe) #create two correlated factors fx <- matrix(c(.9,.8,.7,.85,.75,.65,rep(0,12),.9,.8,.7,.85,.75,.65),ncol=2) Phi <- matrix(c(1,.6,.6,1),2) sim.data <- sim.structure(fx,Phi,n=1000,raw=TRUE) R <- cor(sim.data$observed) Ro <- R[c(1,2,4,5,7,8,10,11),c(1,2,4,5,7,8,10,11)] Roe <- R[c(1,2,4,5,7,8,10,11),c(3,6,9,12)] fo <- fa(Ro,2) fe <- fa.extension(Roe,fo) fa.diagram(fo,fe=fe) #now show how fa.extend works with the same data set #note that we have to make sure that the variables are in the order to do the factor congruence fe2 <- fa.extend(sim.data$observed,2,ov=c(1,2,4,5,7,8,10,11),ev=c(3,6,9,12)) fa.diagram(fe2,main="factor analysis with extension variables") fa2 <- fa(sim.data$observed[,c(1,2,4,5,7,8,10,11,3,6,9,12)],2) factor.congruence(fe2,fa2) summary(fe2) #an example of extending an omega analysis 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(.7,5)) five.factor <- sim.hierarchical(gload,fload,500,TRUE) #create sample data set ss <- c(1,2,3,5,6,7,9,10,11,13,14,15,17,18,19) Ro <- cor(five.factor$observed[,ss]) Re <- cor(five.factor$observed[,ss],five.factor$observed[,-ss]) om5 <-omega(Ro,5) #the omega analysis om.extend <- fa.extension(Re,om5) #the extension analysis om.extend #show it #now, include it in an omega diagram combined.om <- rbind(om5$schmid$sl[,1:ncol(om.extend$loadings)],om.extend$loadings) class(combined.om) <-c("psych","extend") omega.diagram(combined.om,main="Extended Omega") cleanEx() nameEx("fa.lookup") ### * fa.lookup flush(stderr()); flush(stdout()) ### Name: fa.lookup ### Title: A set of functions for factorial and empirical scale ### construction ### Aliases: lookup lookupItems fa.lookup item.lookup keys.lookup ### lookupFromKeys setCorLookup ### Keywords: models multivariate ### ** Examples #Tne following shows how to create a dictionary #first, copy the spreadsheet to the clipboard # bfi.dictionary <- read.clipboard.tab() #read from the clipboard # rownames(bfi.dictionary) <- bfi.dictionary[1] #the first column had the names # bfi.dictionary <- bfi.dictionary[-1] #these are redundant, drop them f5 <- fa(psychTools::bfi,5) m <- colMeans(psychTools::bfi,na.rm=TRUE) item.lookup(f5,m,dictionary=psychTools::bfi.dictionary[2,drop=FALSE]) #just show the item content, not the source of the items fa.lookup(f5,dictionary=psychTools::bfi.dictionary[2]) #show how to use lookupFromKeys bfi.keys <- list(agree=c("-A1","A2","A3","A4","A5"),conscientiousness=c("C1","C2","C3","-C4","-C5"), extraversion=c("-E1","-E2","E3","E4","E5"),neuroticism=c("N1","N2","N3","N4","N5"), openness = c("O1","-O2","O3","O4","-O5")) bfi.over <- scoreOverlap(bfi.keys,bfi) #returns the corrected for overlap values lookupFromKeys(bfi.keys,psychTools::bfi.dictionary,n=5, cors=bfi.over$item.cor) #show the keying information lookupItems("life",psychTools::spi.dictionary) #find those items with "life" in the item 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 fa.parallel.poly plot.poly.parallel ### Keywords: multivariate ### ** Examples #test.data <- Harman74.cor$cov #The 24 variable Holzinger - Harman problem #fa.parallel(test.data,n.obs=145) fa.parallel(Thurstone,n.obs=213) #the 9 variable Thurstone problem #set.seed(123) #minor <- sim.minor(24,4,400) #4 large and 12 minor factors #ffa.parallel(minor$observed) #shows 5 factors and 4 components -- compare with #fa.parallel(minor$observed,SMC=FALSE) #which shows 6 and 4 components factors #a demonstration of parallel analysis of a dichotomous variable #fp <- fa.parallel(psychTools::ability) #use the default Pearson correlation #fpt <- fa.parallel(psychTools::ability,cor="tet") #do a tetrachoric correlation #fpt <- fa.parallel(psychTools::ability,cor="tet",quant=.95) #do a tetrachoric correlation and #use the 95th percentile of the simulated results #apply(fp$values,2,function(x) quantile(x,.95)) #look at the 95th percentile of values #apply(fpt$values,2,function(x) quantile(x,.95)) #look at the 95th percentile of values #describe(fpt$values) #look at all the statistics of the simulated values cleanEx() nameEx("fa.random") ### * fa.random flush(stderr()); flush(stdout()) ### Name: fa.random ### Title: A first approximation to Random Effects Exploratory Factor ### Analysis ### Aliases: fa.random ### Keywords: multivariate models ### ** Examples fa.ab <- fa(psychTools::ability,4,rotate="none") #normal factor analysis fa.ab.ip <- fa.random(psychTools::ability,3,rotate="none") fa.congruence(list(fa.ab,fa.ab.ip,fa.ab.ip$within.r)) cleanEx() nameEx("fa.sort") ### * fa.sort flush(stderr()); flush(stdout()) ### Name: fa.sort ### Title: Sort factor analysis or principal components analysis loadings ### Aliases: fa.sort fa.organize ### Keywords: multivariate ### ** Examples test.simple <- fa(sim.item(16),2) fa.sort(test.simple) fa.organize(test.simple,c(2,1)) #the factors but not the items have been rearranged cleanEx() nameEx("faCor") ### * faCor flush(stderr()); flush(stdout()) ### Name: faCor ### Title: Correlations between two factor analysis solutions ### Aliases: faCor ### Keywords: multivariate models ### ** Examples faCor(Thurstone,nfactors=c(2,3)) #compare two solutions to the Thurstone problem faCor(psychTools::bfi[1:25],nfactors=c(5,5),fm=c("minres","pca")) #compare pca and fa solutions #compare two levels of factor extraction, and then find the correlations of the scores faCor(psychTools::bfi[1:25],nfactors=c(3,5)) #based upon linear algebra f3 <- fa(psychTools::bfi[1:25],3,scores="tenBerge") f5 <- fa(psychTools::bfi[1:25],5 ,scores="tenBerge") cor2(f3$scores,f5$scores) #the correlation between the factor score estimates cleanEx() nameEx("faMulti") ### * faMulti flush(stderr()); flush(stdout()) ### Name: fa.multi ### Title: Multi level (hierarchical) factor analysis ### Aliases: fa.multi fa.multi.diagram ### Keywords: multivariate models ### ** Examples f31 <- fa.multi(Thurstone,3,1) #compare with \code{omega} f31 fa.multi.diagram(f31) cleanEx() nameEx("faRotate") ### * faRotate flush(stderr()); flush(stdout()) ### Name: faRotations ### Title: Multiple rotations of factor loadings to find local minima ### Aliases: faRotations ### Keywords: multivariate models ### ** Examples f5 <- fa(psychTools::bfi[,1:25],5,rotate="none") faRotations(f5,n.rotations=10) #note that the factor analysis needs to not do the rotation faRotations(f5$loadings) #matrix input geo <- faRotations(f5,rotate="geominQ",n.rotation=10) # a popular alternative, but more sensitive to local minima describe(geo$rotation.stats[,1:3]) cleanEx() nameEx("factor.congruence") ### * factor.congruence flush(stderr()); flush(stdout()) ### Name: factor.congruence ### Title: Coefficient of factor congruence ### Aliases: factor.congruence fa.congruence ### Keywords: multivariate models ### ** Examples #factor congruence of factors and components, both rotated #fa <- fa(Harman74.cor$cov,4) #pc <- principal(Harman74.cor$cov,4) #factor.congruence(fa,pc) # RC1 RC3 RC2 RC4 #MR1 0.98 0.41 0.28 0.32 #MR3 0.35 0.96 0.41 0.31 #MR2 0.23 0.16 0.95 0.28 #MR4 0.28 0.38 0.36 0.98 #factor congruence without rotation #fa <- fa(Harman74.cor$cov,4,rotate="none") #pc <- principal(Harman74.cor$cov,4,rotate="none") #factor.congruence(fa,pc) #just show the beween method congruences # PC1 PC2 PC3 PC4 #MR1 1.00 -0.04 -0.06 -0.01 #MR2 0.15 0.97 -0.01 -0.15 #MR3 0.31 0.05 0.94 0.11 #MR4 0.07 0.21 -0.12 0.96 #factor.congruence(list(fa,pc)) #this shows the within method congruence as well # MR1 MR2 MR3 MR4 PC1 PC2 PC3 PC4 #MR1 1.00 0.11 0.25 0.06 1.00 -0.04 -0.06 -0.01 #MR2 0.11 1.00 0.06 0.07 0.15 0.97 -0.01 -0.15 #MR3 0.25 0.06 1.00 0.01 0.31 0.05 0.94 0.11 #MR4 0.06 0.07 0.01 1.00 0.07 0.21 -0.12 0.96 #PC1 1.00 0.15 0.31 0.07 1.00 0.00 0.00 0.00 #PC2 -0.04 0.97 0.05 0.21 0.00 1.00 0.00 0.00 #PC3 -0.06 -0.01 0.94 -0.12 0.00 0.00 1.00 0.00 #PC4 -0.01 -0.15 0.11 0.96 0.00 0.00 0.00 1.00 #pa <- fa(Harman74.cor$cov,4,fm="pa") # 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) # RC1 RC3 RC2 RC4 #MR1 0.99 -0.18 -0.33 -0.34 #MR3 -0.33 0.96 -0.16 -0.43 #MR2 -0.29 -0.46 0.98 -0.21 #MR4 -0.44 -0.30 -0.22 0.98 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) 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 <- fa(Harman74.cor$cov,2,rotate=TRUE) fa2resid <- factor.residuals(Harman74.cor$cov,fa2) fa2resid[1:4,1:4] #residuals with two factors extracted fa4 <- fa(Harman74.cor$cov,4,rotate=TRUE) fa4resid <- factor.residuals(Harman74.cor$cov,fa4) fa4resid[1:4,1:4] #residuals with 4 factors extracted 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 f4<- fa(Harman74.cor$cov,4,rotate="TRUE") f4r45 <- factor.rotate(f4,45,2,3) f4r90 <- factor.rotate(f4r45,45,2,3) print(factor.congruence(f4,f4r45),digits=3) #poor congruence with original print(factor.congruence(f4,f4r90),digits=3) #factor 2 and 3 have been exchanged and 3 flipped #a graphic example data(Harman23.cor) f2 <- fa(Harman23.cor$cov,2,rotate="none") op <- par(mfrow=c(1,2)) cluster.plot(f2,xlim=c(-1,1),ylim=c(-1,1),title="Unrotated ") f2r <- factor.rotate(f2,-33,plot=TRUE,xlim=c(-1,1),ylim=c(-1,1),title="rotated -33 degrees") op <- par(mfrow=c(1,1)) graphics::par(get("par.postscript", pos = 'CheckExEnv')) cleanEx() nameEx("factor.scores") ### * factor.scores flush(stderr()); flush(stdout()) ### Name: factor.scores ### Title: Various ways to estimate factor scores for the factor analysis ### model ### Aliases: factor.scores ### Keywords: multivariate models ### ** Examples f3 <- fa(Thurstone) f3$weights #just the scoring weights f5 <- fa(psychTools::bfi,5) #this does the factor analyis my.scores <- factor.scores(psychTools::bfi,f5, method="tenBerge") #compare the tenBerge factor score correlation to the factor correlations cor(my.scores$scores,use="pairwise") - f5$Phi #compare to the f5$Phi values #compare the default (regression) score correlations to the factor correlations cor(f5$scores,use="pairwise") - f5$Phi #compare to the f5 solution 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 fa.stats ### Keywords: multivariate models ### ** Examples v9 <- sim.hierarchical() f3 <- fa(v9,3) factor.stats(v9,f3,n.obs=500) f3o <- fa(v9,3,fm="pa",rotate="Promax") factor.stats(v9,f3o,n.obs=500) cleanEx() nameEx("factor2cluster") ### * factor2cluster flush(stderr()); flush(stdout()) ### Name: factor2cluster ### Title: Extract cluster definitions from factor loadings ### Aliases: factor2cluster ### Keywords: multivariate models ### ** Examples #matches the factanal example f4 <- fa(Harman74.cor$cov,4,rotate="varimax") factor2cluster(f4) cleanEx() nameEx("fisherz") ### * fisherz flush(stderr()); flush(stdout()) ### Name: fisherz ### Title: Transformations of r, d, and t including Fisher r to z and z to ### r and confidence intervals ### Aliases: fisherz fisherz2r r.con r2c r2t t2r g2r chi2r r2chi cor2cov ### Keywords: multivariate models ### ** Examples n <- 30 r <- seq(0,.9,.1) d <- r2d(r) rc <- matrix(r.con(r,n),ncol=2) t <- r*sqrt(n-2)/sqrt(1-r^2) p <- (1-pt(t,n-2))*2 r1 <- t2r(t,(n-2)) r2 <- d2r(d) chi <- r2chi(r,n) r3 <- chi2r(chi,n) r.rc <- data.frame(r=r,z=fisherz(r),lower=rc[,1],upper=rc[,2],t=t,p=p,d=d, chi2=chi,d2r=r2,t2r=r1,chi2r=r3) round(r.rc,2) cleanEx() nameEx("fparse") ### * fparse flush(stderr()); flush(stdout()) ### Name: fparse ### Title: Parse and exten formula input from a model and return the DV, ### IV, and associated terms. ### Aliases: fparse ### Keywords: utilities ### ** Examples fparse(DV ~ IV1 + IV2 * IV2*IV3 + (IV4) + I(IV5^2) ) #somewhat more complicated fparse(DV1 + DV2 ~ IV1 + IV2 + IV3*IV4 + I(IV5^2) + I(Iv6^2) + (IV7) + (IV8) - IV9) 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 x2[2] <- NA X <- data.frame(x,x2) geometric.mean(x) geometric.mean(x2) geometric.mean(X) geometric.mean(X,na.rm=FALSE) cleanEx() nameEx("glb.algebraic") ### * glb.algebraic flush(stderr()); flush(stdout()) ### Name: glb.algebraic ### Title: Find the greatest lower bound to reliability. ### Aliases: glb.algebraic ### Keywords: multivariate ### ** Examples Cv<-matrix(c(215, 64, 33, 22, 64, 97, 57, 25, 33, 57,103, 36, 22, 25, 36, 77),ncol=4) Cv # covariance matrix of a test with 4 subtests Cr<-cov2cor(Cv) # Correlation matrix of tests if(!require(Rcsdp)) {print("Rcsdp must be installed to find the glb.algebraic")} else { glb.algebraic(Cv) # glb of total score glb.algebraic(Cr) # glb of sum of standardized scores w<-c(1,2,2,1) # glb of weighted total score glb.algebraic(diag(w) %*% Cv %*% diag(w)) alphas <- c(0.8,0,0,0) # Internal consistency of first test is known glb.algebraic(Cv,LoBounds=alphas*diag(Cv)) # Fix all diagonal elements to 1 but the first: lb <- glb.algebraic(Cr,LoBounds=c(0,1,1,1),UpBounds=c(1,1,1,1)) lb$solution[1] # should be the same as the squared mult. corr. smc(Cr)[1] } cleanEx() nameEx("guttman") ### * guttman flush(stderr()); flush(stdout()) ### Name: splitHalf ### Title: Alternative estimates of test reliabiity ### Aliases: splitHalf guttman tenberge glb glb.fa ### Keywords: multivariate ### ** Examples data(attitude) splitHalf(attitude) splitHalf(attitude,covar=TRUE) #do it on the covariances temp <- splitHalf(attitude,raw=TRUE) temp$ci #to show the confidence intervals, you need to specify that raw=TRUE glb(attitude) glb.fa(attitude) if(require(Rcsdp)) {glb.algebraic(cor(attitude)) } guttman(attitude) #to show the histogram of all possible splits for the ability test #sp <- splitHalf(psychTools::ability,raw=TRUE) #this saves the results #hist(sp$raw,breaks=101,ylab="SplitHalf reliability",main="SplitHalf # reliabilities of a test with 16 ability items") sp <- splitHalf(psychTools::bfi[1:10],key=c(1,9,10)) 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 x2[2] <- NA y <- x - 1 X <- data.frame(x,x2,y) harmonic.mean(x) harmonic.mean(x2) harmonic.mean(X) harmonic.mean(X,na.rm=FALSE) harmonic.mean(X,zero=FALSE) cleanEx() nameEx("headtail") ### * headtail flush(stderr()); flush(stdout()) ### Name: headTail ### Title: Combine calls to head and tail ### Aliases: headtail headTail topBottom quickView ### Keywords: multivariate ### ** Examples headTail(psychTools::iqitems,4,8,to=6) #the first 4 and last 6 items from 1 to 6 topBottom(psychTools::ability,from =2, to = 6) #the first and last 4 items from 2 to 6 headTail(psychTools::bfi,top=4, bottom=4,from =6,to=10) #the first and last 4 from 6 to 10 #not shown #quickView(ability,hlength=10,tlength=10) #this loads a spreadsheet like table cleanEx() nameEx("iclust.diagram") ### * iclust.diagram flush(stderr()); flush(stdout()) ### Name: iclust.diagram ### Title: Draw an ICLUST hierarchical cluster structure diagram ### Aliases: iclust.diagram ICLUST.diagram ### Keywords: multivariate cluster hplot ### ** Examples v9 <- sim.hierarchical() v9c <- ICLUST(v9) test.data <- Harman74.cor$cov ic.out <- ICLUST(test.data) #now show how to relabel clusters ic.bfi <- iclust(psychTools::bfi[1:25],beta=3) #find the clusters cluster.names <- rownames(ic.bfi$results) #get the old names #change the names to the desired ones cluster.names[c(16,19,18,15,20)] <- c("Neuroticism","Extra-Open","Agreeableness", "Conscientiousness","Open") #now show the new names iclust.diagram(ic.bfi,cluster.names=cluster.names,min.size=4,e.size=1.75) 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 interp.median(c(1,2,2,5)) interp.quantiles(c(1,2,2,5),.25) x <- sample(10,100,TRUE) interp.quartiles(x) # 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,4,4,1,2,3,3,3,3,4,4,4,4,5,1,5,3,3,3,3,4,4,4,4,4) 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(psychTools::galton) galton <- psychTools::galton interp.median(galton) interp.qplot.by(galton$child,galton$parent,ylab="child height" ,xlab="Mid parent height") cleanEx() nameEx("irt.fa") ### * irt.fa flush(stderr()); flush(stdout()) ### Name: irt.fa ### Title: Item Response Analysis by Exploratory Factor Analysis of ### tetrachoric/polychoric correlations ### Aliases: irt.fa irt.select fa2irt ### Keywords: multivariate models ### ** Examples ## Not run: ##D set.seed(17) ##D d9 <- sim.irt(9,1000,-2.5,2.5,mod="normal") #dichotomous items ##D test <- irt.fa(d9$items) ##D test ##D op <- par(mfrow=c(3,1)) ##D plot(test,type="ICC") ##D plot(test,type="IIC") ##D plot(test,type="test") ##D par(op) ##D set.seed(17) ##D items <- sim.congeneric(N=500,short=FALSE,categorical=TRUE) #500 responses to 4 discrete items ##D d4 <- irt.fa(items$observed) #item response analysis of congeneric measures ##D d4 #show just the irt output ##D d4$fa #show just the factor analysis output ##D ##D ##D op <- par(mfrow=c(2,2)) ##D plot(d4,type="ICC") ##D par(op) ##D ##D ##D #using the iq data set for an example of real items ##D #first need to convert the responses to tf ##D data(iqitems) ##D iq.keys <- c(4,4,4, 6, 6,3,4,4, 5,2,2,4, 3,2,6,7) ##D ##D iq.tf <- score.multiple.choice(iq.keys,psychTools::iqitems,score=FALSE) #just the responses ##D iq.irt <- irt.fa(iq.tf) ##D print(iq.irt,short=FALSE) #show the IRT as well as factor analysis output ##D p.iq <- plot(iq.irt) #save the invisible summary table ##D p.iq #show the summary table of information by ability level ##D #select a subset of these variables ##D small.iq.irt <- irt.select(iq.irt,c(1,5,9,10,11,13)) ##D small.irt <- irt.fa(small.iq.irt) ##D plot(small.irt) ##D #find the information for three subset of iq items ##D keys <- make.keys(16,list(all=1:16,some=c(1,5,9,10,11,13),others=c(1:5))) ##D plot(iq.irt,keys=keys) ## End(Not run) #compare output to the ltm package or Kamata and Bauer -- these are in logistic units ls <- irt.fa(lsat6) #library(ltm) # lsat.ltm <- ltm(lsat6~z1) # round(coefficients(lsat.ltm)/1.702,3) #convert to normal (approximation) # # Dffclt Dscrmn #Q1 -1.974 0.485 #Q2 -0.805 0.425 #Q3 -0.164 0.523 #Q4 -1.096 0.405 #Q5 -1.835 0.386 #Normal results ("Standardized and Marginal")(from Akihito Kamata ) #Item discrim tau # 1 0.4169 -1.5520 # 2 0.4333 -0.5999 # 3 0.5373 -0.1512 # 4 0.4044 -0.7723 # 5 0.3587 -1.1966 #compare to ls #Normal results ("Standardized and conditional") (from Akihito Kamata ) #item discrim tau # 1 0.3848 -1.4325 # 2 0.3976 -0.5505 # 3 0.4733 -0.1332 # 4 0.3749 -0.7159 # 5 0.3377 -1.1264 #compare to ls$fa and ls$tau #Kamata and Bauer (2008) logistic estimates #1 0.826 2.773 #2 0.723 0.990 #3 0.891 0.249 #4 0.688 1.285 #5 0.657 2.053 cleanEx() nameEx("irt.responses") ### * irt.responses flush(stderr()); flush(stdout()) ### Name: irt.responses ### Title: Plot probability of multiple choice responses as a function of a ### latent trait ### Aliases: irt.responses ### Keywords: multivariate models ### ** Examples data(psychTools::iqitems) iq.keys <- c(4,4,4, 6,6,3,4,4, 5,2,2,4, 3,2,6,7) scores <- score.multiple.choice(iq.keys,psychTools::iqitems,score=TRUE,short=FALSE) #note that for speed we can just do this on simple item counts rather # than IRT based scores. op <- par(mfrow=c(2,2)) #set this to see the output for multiple items irt.responses(scores$scores,psychTools::iqitems[1:4],breaks=11) op <- par(op) graphics::par(get("par.postscript", pos = 'CheckExEnv')) cleanEx() nameEx("kaiser") ### * kaiser flush(stderr()); flush(stdout()) ### Name: kaiser ### Title: Apply the Kaiser normalization when rotating factors ### Aliases: kaiser ### Keywords: multivariate models ### ** Examples f3 <- fa(Thurstone,3) f3n <- kaiser(fa(Thurstone,3,rotate="none")) f3p <- kaiser(fa(Thurstone,3,rotate="none"),rotate="Promax",m=3) factor.congruence(list(f3,f3n,f3p)) cleanEx() nameEx("kappa") ### * kappa flush(stderr()); flush(stdout()) ### Name: cohen.kappa ### Title: Find Cohen's kappa and weighted kappa coefficients for ### correlation of two raters ### Aliases: wkappa cohen.kappa ### Keywords: multivariate ### ** Examples #rating data (with thanks to Tim Bates) rater1 = c(1,2,3,4,5,6,7,8,9) # rater one's ratings rater2 = c(1,3,1,6,1,5,5,6,7) # rater one's ratings cohen.kappa(x=cbind(rater1,rater2)) #data matrix taken from Cohen cohen <- matrix(c( 0.44, 0.07, 0.09, 0.05, 0.20, 0.05, 0.01, 0.03, 0.06),ncol=3,byrow=TRUE) #cohen.weights weight differences cohen.weights <- matrix(c( 0,1,3, 1,0,6, 3,6,0),ncol=3) cohen.kappa(cohen,cohen.weights,n.obs=200) #cohen reports .492 and .348 #another set of weights #what if the weights are non-symmetric wc <- matrix(c( 0,1,4, 1,0,6, 2,2,0),ncol=3,byrow=TRUE) cohen.kappa(cohen,wc) #Cohen reports kw = .353 cohen.kappa(cohen,n.obs=200) #this uses the squared weights fleiss.cohen <- 1 - cohen.weights/9 cohen.kappa(cohen,fleiss.cohen,n.obs=200) #however, Fleiss, Cohen and Everitt weight similarities fleiss <- matrix(c( 106, 10,4, 22,28, 10, 2, 12, 6),ncol=3,byrow=TRUE) #Fleiss weights the similarities weights <- matrix(c( 1.0000, 0.0000, 0.4444, 0.0000, 1.0000, 0.6667, 0.4444, 0.6667, 1.0000),ncol=3) cohen.kappa(fleiss,weights,n.obs=200) #another example is comparing the scores of two sets of twins #data may be a 2 column matrix #compare weighted and unweighted #also look at the ICC for this data set. twins <- matrix(c( 1, 2, 2, 3, 3, 4, 5, 6, 6, 7), ncol=2,byrow=TRUE) cohen.kappa(twins) #data may be explicitly categorical x <- c("red","yellow","blue","red") y <- c("red", "blue", "blue" ,"red") xy.df <- data.frame(x,y) ck <- cohen.kappa(xy.df) ck ck$agree #Example for specifying levels #The problem of missing categories (from Amy Finnegan) #We need to specify all the categories possible using the levels option numbers <- data.frame(rater1=c(6,3,7,8,7), rater2=c(6,1,8,5,10)) cohen.kappa(numbers) #compare with the next analysis cohen.kappa(numbers,levels=1:10) #specify the number of levels # these leads to slightly higher weighted kappa #finally, input can be a data.frame of ratings from more than two raters ratings <- matrix(rep(1:5,4),ncol=4) ratings[1,2] <- ratings[2,3] <- ratings[3,4] <- NA ratings[2,1] <- ratings[3,2] <- ratings[4,3] <- 1 ck <- cohen.kappa(ratings) ck #just show the raw and weighted kappas print(ck, all=TRUE) #show the confidence intervals as well #In the case of confidence intervals being artificially truncated to +/- 1, it is #helpful to compare the results of a boot strap resample #ck.boot <-function(x,s=1:nrow(x)) {cohen.kappa(x[s,])$kappa} #library(boot) #ckb <- boot(x,ck.boot,R=1000) #hist(ckb$t) cleanEx() nameEx("logistic") ### * logistic flush(stderr()); flush(stdout()) ### Name: logistic ### Title: Logistic transform from x to p and logit transform from p to x ### Aliases: logistic logit logistic.grm ### Keywords: multivariate ### ** Examples curve(logistic(x,a=1.702),-3,3,ylab="Probability of x", main="Logistic transform of x",xlab="z score units") #logistic with a=1.702 is almost the same as pnorm curve(pnorm(x),add=TRUE,lty="dashed") curve(logistic(x),add=TRUE) text(2,.8, expression(alpha ==1)) text(2,1.0,expression(alpha==1.7)) curve(logistic(x),-4,4,ylab="Probability of x", main = "Logistic transform of x in logit units",xlab="logits") curve(logistic(x,d=-1),add=TRUE) curve(logistic(x,d=1),add=TRUE) curve(logistic(x,c=.2),add=TRUE,lty="dashed") text(1.3,.5,"d=1") text(.3,.5,"d=0") text(-1.5,.5,"d=-1") text(-3,.3,"c=.2") #demo of graded response model curve(logistic.grm(x,r=1),-4,4,ylim=c(0,1),main="Five level response scale", ylab="Probability of endorsement",xlab="Latent attribute on logit scale") curve(logistic.grm(x,r=2),add=TRUE) curve(logistic.grm(x,r=3),add=TRUE) curve(logistic.grm(x,r=4),add=TRUE) curve(logistic.grm(x,r=5),add=TRUE) text(-2.,.5,1) text(-1.,.4,2) text(0,.4,3) text(1.,.4,4) text(2.,.4,5) cleanEx() nameEx("lowerUpper") ### * lowerUpper flush(stderr()); flush(stdout()) ### Name: lowerUpper ### Title: Combine two square matrices to have a lower off diagonal for ### one, upper off diagonal for the other ### Aliases: lowerUpper ### Keywords: multivariate ### ** Examples b1 <- Bechtoldt.1 b2 <- Bechtoldt.2 b12 <- lowerUpper(b1,b2) cor.plot(b12) diff12 <- lowerUpper(b1,b2,diff=TRUE) corPlot(t(diff12),numbers=TRUE,main="Bechtoldt1 and the differences from Bechtoldt2") #Compare r and partial r lower <- lowerCor(sat.act) upper <- partial.r(sat.act) both = lowerUpper(lower,upper) corPlot(both,numbers=TRUE,main="r and partial r for the sat.act data set") #now show the differences both = lowerUpper(lower,upper,diff=TRUE) corPlot(both,numbers=TRUE,main="Differences between r and partial r for the sat.act data set") 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 keys2list selectFromKeys makePositiveKeys ### Keywords: multivariate models ### ** Examples data(attitude) #specify the items by location 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 #now, undo this new.keys.list <- keys2list(keys) #note, these are now given the variable names select <- selectFromKeys(key.list) #scores <- score.items(keys,attitude) #scores # data(psychTools::bfi) #first create the keys by location (the conventional way) 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(psychTools::bfi)[1:25]) new.keys.list <- keys2list(keys) #these will be in the form of variable names #alternatively, create by a mixture of names and locations keys.list <- list(agree=c("-A1","A2","A3","A4","A5"), conscientious=c("C1","C2","C2","-C4","-C5"),extraversion=c("-E1","-E2","E3","E4","E5"), neuroticism=c(16:20),openness = c(21,-22,23,24,-25)) keys <- make.keys(psychTools::bfi, keys.list) #specify the data file to be scored (bfi) #or keys <- make.keys(colnames(psychTools::bfi),keys.list) #specify the names of the variables #to be used #or #specify the number of variables to be scored and their names in all cases keys <- make.keys(28,keys.list,colnames(psychTools::bfi)) scores <- scoreItems(keys,psychTools::bfi) summary(scores) cleanEx() nameEx("manhattan") ### * manhattan flush(stderr()); flush(stdout()) ### Name: manhattan ### Title: "Manhattan" plots of correlations with a set of criteria. ### Aliases: manhattan ### Keywords: multivariate hplot ### ** Examples op <- par(mfrow=(c(2,3))) #we want to compare two different sets of plots manhattan(psychTools::bfi[1:25],psychTools::bfi[26:28] ,labels=colnames(psychTools::bfi)[1:25], dictionary=psychTools::bfi.dictionary) manhattan(psychTools::bfi[1:25],psychTools::bfi[26:28],log.p=TRUE, dictionary=psychTools::bfi.dictionary) #Do it again, but now show items by the keys.list bfi.keys <- list(agree=c("-A1","A2","A3","A4","A5"),conscientious=c("C1","C2","C3","-C4","-C5"), extraversion=c("-E1","-E2","E3","E4","E5"),neuroticism=c("N1","N2","N3","N4","N5"), openness = c("O1","-O2","O3","O4","-O5")) man <- manhattan(psychTools::bfi[1:25],psychTools::bfi[26:28],keys=bfi.keys, dictionary=psychTools::bfi.dictionary[1:2]) manhattan(psychTools::bfi[1:25],psychTools::bfi[26:28],keys=bfi.keys,log.p=TRUE, dictionary=psychTools::bfi.dictionary[1:2]) #Alternatively, use a matrix as input R <-cor(psychTools::bfi[1:25],psychTools::bfi[26:28],use="pairwise") manhattan(R,cs(gender,education,age),keys=bfi.keys, dictionary=psychTools::bfi.dictionary[1:2], raw=FALSE,abs=FALSE) par <- op psychTools::dfOrder(man,1,ascending=FALSE) #print out the items sorted on gender graphics::par(get("par.postscript", pos = 'CheckExEnv')) 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 matSort ### Keywords: multivariate models ### ** Examples data(Bechtoldt.1) sorted <- mat.sort(Bechtoldt.1,fa(Bechtoldt.1,5)) corPlot(sorted,xlas=2) #vertical xaxis names 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 z #compare with outer(x,-x,FUN="+") x <- matrix(seq(1,6),ncol=2) y <- matrix(seq(1,10),nrow=2) z <- x %+% y x y z #but compare this with outer(x ,y,FUN="+") cleanEx() nameEx("mediate") ### * mediate flush(stderr()); flush(stdout()) ### Name: mediate ### Title: Estimate and display direct and indirect effects of mediators ### and moderator in path models ### Aliases: mediate mediate.diagram moderate.diagram ### Keywords: multivariate models ### ** Examples # A simple mediation example is the Tal_Or data set (pmi for Hayes) #The pmi data set from Hayes is available as the Tal_Or data set. mod4 <- mediate(reaction ~ cond + (pmi), data =Tal_Or,n.iter=50) summary(mod4) #Two mediators (from Hayes model 6 (chapter 5)) mod6 <- mediate(reaction ~ cond + (pmi) + (import), data =Tal_Or,n.iter=50) summary(mod6) #Moderated mediation is done for the Garcia (Garcia, 2010) data set. # (see Hayes, 2013 for the protest data set #n.iter set to 50 (instead of default of 5000) for speed of example #no mediation, just an interaction mod7 <- mediate(liking ~ sexism * prot2 , data=Garcia, n.iter = 50) summary(mod7) data(GSBE) #The Garcia et al data set (aka GSBE) mod11.4 <- mediate(liking ~ sexism * prot2 + (respappr), data=Garcia, n.iter = 50,zero=FALSE) #to match Hayes summary(mod11.4) #to see this interaction graphically, run the examples in ?Garcia #data from Preacher and Hayes (2004) sobel <- structure(list(SATIS = c(-0.59, 1.3, 0.02, 0.01, 0.79, -0.35, -0.03, 1.75, -0.8, -1.2, -1.27, 0.7, -1.59, 0.68, -0.39, 1.33, -1.59, 1.34, 0.1, 0.05, 0.66, 0.56, 0.85, 0.88, 0.14, -0.72, 0.84, -1.13, -0.13, 0.2), THERAPY = structure(c(0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0), value.labels = structure(c(1, 0), .Names = c("cognitive", "standard"))), ATTRIB = c(-1.17, 0.04, 0.58, -0.23, 0.62, -0.26, -0.28, 0.52, 0.34, -0.09, -1.09, 1.05, -1.84, -0.95, 0.15, 0.07, -0.1, 2.35, 0.75, 0.49, 0.67, 1.21, 0.31, 1.97, -0.94, 0.11, -0.54, -0.23, 0.05, -1.07)), .Names = c("SATIS", "THERAPY", "ATTRIB" ), row.names = c(NA, -30L), class = "data.frame", variable.labels = structure(c("Satisfaction", "Therapy", "Attributional Positivity"), .Names = c("SATIS", "THERAPY", "ATTRIB"))) #n.iter set to 50 (instead of default of 5000) for speed of example #There are several forms of input. The original specified y, x , and the mediator #mediate(1,2,3,sobel,n.iter=50) #The example in Preacher and Hayes #As of October, 2017 we can specify this in a formula mode mediate (SATIS ~ THERAPY + (ATTRIB),data=sobel, n.iter=50) #specify the mediator by # adding parentheses #A.C. Kerchoff, (1974) Ambition and Attainment: A Study of Four Samples of American Boys. #Data from sem package taken from Kerckhoff (and in turn, from Lisrel manual) R.kerch <- structure(list(Intelligence = c(1, -0.1, 0.277, 0.25, 0.572, 0.489, 0.335), Siblings = c(-0.1, 1, -0.152, -0.108, -0.105, -0.213, -0.153), FatherEd = c(0.277, -0.152, 1, 0.611, 0.294, 0.446, 0.303), FatherOcc = c(0.25, -0.108, 0.611, 1, 0.248, 0.41, 0.331), Grades = c(0.572, -0.105, 0.294, 0.248, 1, 0.597, 0.478 ), EducExp = c(0.489, -0.213, 0.446, 0.41, 0.597, 1, 0.651), OccupAsp = c(0.335, -0.153, 0.303, 0.331, 0.478, 0.651, 1 )), .Names = c("Intelligence", "Siblings", "FatherEd", "FatherOcc", "Grades", "EducExp", "OccupAsp"), class = "data.frame", row.names = c("Intelligence", "Siblings", "FatherEd", "FatherOcc", "Grades", "EducExp", "OccupAsp" )) #n.iter set to 50 (instead of default of 5000) for speed of demo #mod.k <- mediate("OccupAsp","Intelligence",m= c(2:5),data=R.kerch,n.obs=767,n.iter=50) #new style mod.k <- mediate(OccupAsp ~ Intelligence + (Siblings) + (FatherEd) + (FatherOcc) + (Grades), data = R.kerch, n.obs=767, n.iter=50) mediate.diagram(mod.k) #print the path values mod.k #Compare the following solution to the path coefficients found by the sem package #mod.k2 <- mediate(y="OccupAsp",x=c("Intelligence","Siblings","FatherEd","FatherOcc"), # m= c(5:6),data=R.kerch,n.obs=767,n.iter=50) #new format mod.k2 <- mediate(OccupAsp ~ Intelligence + Siblings + FatherEd + FatherOcc + (Grades) + (EducExp),data=R.kerch, n.obs=767, n.iter=50) mediate.diagram(mod.k2,show.c=FALSE) #simpler output #print the path values mod.k2 #Several interesting test cases are taken from analyses of the Spengler data set #This is temporarily added to psych from psychTools to help build for CRAN #Although the sample sizes are actually very large in the first wave, I use the #sample sizes from the last wave #We set the n.iter to be 50 instead of the default value of 5,000 if(require("psychTools")) { mod1 <- mediate(Income.50 ~ IQ + Parental+ (Ed.11) ,data=Spengler, n.obs = 1952, n.iter=50) mod2 <- mediate(Income.50 ~ IQ + Parental+ (Ed.11) + (Income.11) ,data=Spengler,n.obs = 1952, n.iter=50) #Now, compare these models anova(mod1,mod2) #Current version does not support two DVs #mod22 <- mediate(Income.50 + Educ.50 ~ IQ + Parental+ (Ed.11) + (Income.11) # ,data=Spengler,n.obs = 1952, n.iter=50) } cleanEx() nameEx("misc") ### * misc flush(stderr()); flush(stdout()) ### Name: psych.misc ### Title: Miscellaneous helper functions for the psych package ### Aliases: psych.misc misc tableF lowerCor lowerMat progressBar reflect ### shannon test.all cor2 levels2numeric char2numeric nchar2numeric ### isCorrelation isCovariance fromTo cs acs ### Keywords: multivariate ### ** Examples lowerMat(Thurstone) lb <- lowerCor(psychTools::bfi[1:10]) #finds and prints the lower correlation matrix, # returns the square matrix. #fiml <- corFiml(psychTools::bfi[1:10]) #FIML correlations require lavaan package #lowerMat(fiml) #to get pretty output f3 <- fa(Thurstone,3) f3r <- reflect(f3,2) #reflect the second factor #find the complexity of the response patterns of the iqitems. round(shannon(psychTools::iqitems),2) #test.all('BinNor') #Does the BinNor package work when we are using other packages bestItems(lb,"A3",cut=.1,dictionary=psychTools::bfi.dictionary[1:2],raw=FALSE) #to make this a latex table #df2latex(bestItems(lb,2,cut=.2)) # data(psychTools::bfi.dictionary) f2 <- fa(psychTools::bfi[1:10],2) fa.lookup(f2,psychTools::bfi.dictionary) sa1 <-sat.act[1:2] sa2 <- sat.act[3:4] sa3 <- sat.act[5:6] cor2(sa1,sa2) cor2(list(sa1,sa2)) #show within set and between set cors cor2(list(sa1,sa2,sa3)) lowerCor(fromTo(sat.act,"ACT","SATQ")) #show some correlations vect <- cs(ACT,SATQ) #skip the quotes vect #they are in this vector #to combine longer terms vect <- cs("Here is a longish",vector, that, we ,"want to combine", into, several) vect temp <- acs("Here is a longish",vector, that, we ,"want to combine", into, one) temp lowerCor(fromTo(sat.act,cs(ACT,SATQ))) cleanEx() nameEx("mixed.cor") ### * mixed.cor flush(stderr()); flush(stdout()) ### Name: mixedCor ### Title: Find correlations for mixtures of continuous, polytomous, and ### dichotomous variables ### Aliases: mixedCor mixed.cor ### Keywords: multivariate models ### ** Examples data(bfi) r <- mixedCor(data=psychTools::bfi[,c(1:5,26,28)]) r #this is the same as r <- mixedCor(data=psychTools::bfi,p=1:5,c=28,d=26) r #note how the variable order reflects the original order in the data #compare to raw Pearson #note that the biserials and polychorics are not attenuated rp <- cor(psychTools::bfi[c(1:5,26,28)],use="pairwise") lowerMat(rp) cleanEx() nameEx("mssd") ### * mssd flush(stderr()); flush(stdout()) ### Name: mssd ### Title: Find von Neuman's Mean Square of Successive Differences ### Aliases: mssd rmssd autoR ### Keywords: multivariate models ### ** Examples t <- seq(-pi, pi, .1) trial <- 1: length(t) gr <- trial %% 8 c <- cos(t) ts <- sample(t,length(t)) cs <- cos(ts) x.df <- data.frame(trial,gr,t,c,ts,cs) rmssd(x.df) rmssd(x.df,gr) autoR(x.df,gr) describe(x.df) #pairs.panels(x.df) #mlPlot(x.df,grp="gr",Time="t",items=c(4:6)) 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 histBy ### Keywords: multivariate hplot ### ** Examples multi.hist(sat.act) multi.hist(sat.act,bcol="red") multi.hist(sat.act,dcol="blue") #make both lines blue multi.hist(sat.act,dcol= c("blue","red"),dlty=c("dotted", "solid")) multi.hist(sat.act,freq=TRUE) #show the frequency plot multi.hist(sat.act,nrow=2) histBy(sat.act,"SATQ","gender") #input by variable names histBy(SATQ~ gender, data=sat.act) #formula input cleanEx() nameEx("multilevel.reliability") ### * multilevel.reliability flush(stderr()); flush(stdout()) ### Name: multilevel.reliability ### Title: Find and plot various reliability/gneralizability coefficients ### for multilevel data ### Aliases: mlr multilevel.reliability mlArrange mlPlot ### Keywords: multivariate models ### ** Examples #data from Shrout and Lane, 2012. shrout <- structure(list(Person = c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), Time = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L), Item1 = c(2L, 3L, 6L, 3L, 7L, 3L, 5L, 6L, 3L, 8L, 4L, 4L, 7L, 5L, 6L, 1L, 5L, 8L, 8L, 6L), Item2 = c(3L, 4L, 6L, 4L, 8L, 3L, 7L, 7L, 5L, 8L, 2L, 6L, 8L, 6L, 7L, 3L, 9L, 9L, 7L, 8L ), Item3 = c(6L, 4L, 5L, 3L, 7L, 4L, 7L, 8L, 9L, 9L, 5L, 7L, 9L, 7L, 8L, 4L, 7L, 9L, 9L, 6L)), .Names = c("Person", "Time", "Item1", "Item2", "Item3"), class = "data.frame", row.names = c(NA, -20L)) #make shrout super wide #Xwide <- reshape(shrout,v.names=c("Item1","Item2","Item3"),timevar="Time", #direction="wide",idvar="Person") #add more helpful Names #colnames(Xwide ) <- c("Person",c(paste0("Item",1:3,".T",1),paste0("Item",1:3,".T",2), #paste0("Item",1:3,".T",3),paste0("Item",1:3,".T",4))) #make superwide into normal form (i.e., just return it to the original shrout data #Xlong <-Xlong <- reshape(Xwide,idvar="Person",2:13) #Now use these data for a multilevel repliability study, use the normal wide form output mg <- mlr(shrout,grp="Person",Time="Time",items=3:5) #which is the same as #mg <- multilevel.reliability(shrout,grp="Person",Time="Time",items= # c("Item1","Item2","Item3"),plot=TRUE) #to show the lattice plot by subjects, set plot = TRUE #Alternatively for long input (returned in this case from the prior run) mlr(mg$long,grp="id",Time ="time",items="items", values="values",long=TRUE) #example of mlArrange #First, add two new columns to shrout and #then convert to long output using mlArrange total <- rowSums(shrout[3:5]) caseid <- rep(paste0("ID",1:5),4) new.shrout <- cbind(shrout,total=total,case=caseid) #now convert to long new.long <- mlArrange(new.shrout,grp="Person",Time="Time",items =3:5,extra=6:7) headTail(new.long,6,6) cleanEx() nameEx("omega") ### * omega flush(stderr()); flush(stdout()) ### Name: omega ### Title: Calculate McDonald's omega estimates of general and total factor ### saturation ### Aliases: omega omegaSem omegaFromSem omegah omegaDirect directSl ### Keywords: multivariate models ### ** Examples ## Not run: ##D test.data <- Harman74.cor$cov ##D # if(!require(GPArotation)) {message("Omega requires GPA rotation" )} else { ##D my.omega <- omega(test.data) ##D print(my.omega,digits=2) ##D #} ##D ##D #create 9 variables with a hierarchical structure ##D v9 <- sim.hierarchical() ##D #with correlations of ##D round(v9,2) ##D #find omega ##D v9.omega <- omega(v9,digits=2) ##D v9.omega ##D ##D #create 8 items with a two factor solution, showing the use of the flip option ##D sim2 <- item.sim(8) ##D omega(sim2) #an example of misidentification-- remember to look at the loadings matrices. ##D omega(sim2,2) #this shows that in fact there is no general factor ##D omega(sim2,2,option="first") #but, if we define one of the two group factors ##D #as a general factor, we get a falsely high omega ##D #apply omega to analyze 6 mental ability tests ##D data(ability.cov) #has a covariance matrix ##D omega(ability.cov$cov) ##D ##D #om <- omega(Thurstone) ##D #round(om$omega.group,2) ##D #round(om$omega.group[2]/om$omega.group[1],2) #fraction of reliable that is general variance ##D # round(om$omega.group[3]/om$omega.group[1],2) #fraction of reliable that is group variance ##D ##D #To find factor score estimates for the hierarchical model it is necessary to ##D #do two extra steps. ##D ##D #Consider the case of the raw data in an object data. (An example from simulation) ##D # set.seed(42) ##D # gload <- matrix(c(.9,.8,.7),nrow=3) ##D # fload <- matrix(c(.8,.7,.6,rep(0,9),.7,.6,.5,rep(0,9),.7,.6,.4), ncol=3) ##D # data <- sim.hierarchical(gload=gload,fload=fload, n=100000, raw=TRUE) ##D # ##D # f3 <- fa(data$observed,3,scores="tenBerge", oblique.scores=TRUE) ##D # f1 <- fa(f3$scores) ##D ##D # om <- omega(data$observed,sl=FALSE) #draw the hierarchical figure ##D # The scores from om are based upon the Schmid-Leiman factors and although the g factor ##D # is identical, the group factors are not. ##D # This is seen in the following correlation matrix ##D # hier.scores <- cbind(om$scores,f1$scores,f3$scores) ##D # lowerCor(hier.scores) ##D # ##D #this next set of examples require lavaan ##D #jensen <- sim.hierarchical() #create a hierarchical structure (same as v9 above) ##D #om.jen <- omegaSem(jensen,lavaan=TRUE) #do the exploratory omega with confirmatory as well ##D #lav.mod <- om.jen$omegaSem$model$lavaan #get the lavaan code or create it yourself ##D # lav.mod <- 'g =~ +V1+V2+V3+V4+V5+V6+V7+V8+V9 ##D # F1=~ + V1 + V2 + V3 ##D # F2=~ + V4 + V5 + V6 ##D # F3=~ + V7 + V8 + V9 ' ##D #lav.jen <- cfa(lav.mod,sample.cov=jensen,sample.nobs=500,orthogonal=TRUE,std.lv=TRUE) ##D # omegaFromSem(lav.jen,jensen) ##D #the directSl solution ##D #direct.jen <- directSl(jen) ##D #direct.jen ##D ##D #try a one factor solution -- this is not recommended, but sometimes done ##D #it will just give omega_total ##D # lav.mod.1 <- 'g =~ +V1+V2+V3+V4+V5+V6+V7+V8+V9 ' ##D #lav.jen.1<- cfa(lav.mod.1,sample.cov=jensen,sample.nobs=500,orthogonal=TRUE,std.lv=TRUE) ##D # omegaFromSem(lav.jen.1,jensen) ##D ##D ##D ## End(Not run) cleanEx() nameEx("omega.graph") ### * omega.graph flush(stderr()); flush(stdout()) ### Name: omega.graph ### Title: Graph hierarchical factor structures ### Aliases: omega.diagram omega.graph ### Keywords: multivariate ### ** Examples #24 mental tests from Holzinger-Swineford-Harman if(require(GPArotation) ) {om24 <- omega(Harman74.cor$cov,4) } #run omega # #example hierarchical structure from Jensen and Weng if(require(GPArotation) ) {jen.omega <- omega(make.hierarchical())} cleanEx() nameEx("outlier") ### * outlier flush(stderr()); flush(stdout()) ### Name: outlier ### Title: Find and graph Mahalanobis squared distances to detect outliers ### Aliases: outlier ### Keywords: multivariate models ### ** Examples #first, just find and graph the outliers d2 <- outlier(sat.act) #combine with the data frame and plot it with the outliers highlighted in blue sat.d2 <- data.frame(sat.act,d2) pairs.panels(sat.d2,bg=c("yellow","blue")[(d2 > 25)+1],pch=21) 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 p.rep.f(9.0,98) #probability of replicating a result with F = 9.0 with 98 df p.rep.r(.4,50) #probability of replicating a result if r =.4 with n = 50 p.rep.t(3,98) #probability of replicating a result if t = 3 with df =98 p.rep.t(2.14,84,14) #effect of equal sample sizes (see Rosnow et al.) 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 paired.r(.5,.3,NULL,100) #independent correlations same sample size paired.r(.5,.3,NULL, 100, 64) # independent correlations, different sample sizes 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 panel.smoother ### Keywords: multivariate hplot ### ** Examples pairs.panels(attitude) #see the graphics window data(iris) pairs.panels(iris[1:4],bg=c("red","yellow","blue")[iris$Species], pch=21,main="Fisher Iris data by Species") #to show color grouping pairs.panels(iris[1:4],bg=c("red","yellow","blue")[iris$Species], pch=21+as.numeric(iris$Species),main="Fisher Iris data by Species",hist.col="red") #to show changing the diagonal #to show 'significance' pairs.panels(iris[1:4],bg=c("red","yellow","blue")[iris$Species], pch=21+as.numeric(iris$Species),main="Fisher Iris data by Species",hist.col="red",stars=TRUE) #demonstrate not showing the data points data(sat.act) pairs.panels(sat.act,show.points=FALSE) #better yet is to show the points as a period pairs.panels(sat.act,pch=".") #show many variables with 0 gap between scatterplots # data(bfi) # pairs.panels(psychTools::bfi,show.points=FALSE,gap=0) #plot raw data points and then the weighted correlations. #output from statsBy sb <- statsBy(sat.act,"education") pairs.panels(sb$mean,wt=sb$n) #report the weighted correlations #compare with pairs.panels(sb$mean) #unweighed correlations cleanEx() nameEx("parcels") ### * parcels flush(stderr()); flush(stdout()) ### Name: parcels ### Title: Find miniscales (parcels) of size 2 or 3 from a set of items ### Aliases: parcels keysort ### Keywords: multivariate ### ** Examples parcels(Thurstone) keys <- parcels(psychTools::bfi) keys <- keysort(keys) score.items(keys,psychTools::bfi) 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 lowerMat(jen[1:5,1:5]) par.r <- partial.r(jen,c(1,3,5),c(2,4)) lowerMat(par.r) #or R <- jen[1:5,1:5] par.r <- partial.r(R, y = cs(V2,V4)) lowerMat(par.r) cp <- corr.p(par.r,n=98) #assumes the jen data based upon n =100. print(cp,short=FALSE) #show the confidence intervals as well #partial all from all correlations. lowerMat(partial.r(jen)) #Consider the Tal.Or data set. lowerCor(Tal.Or) #partial gender and age from these relations (they hardly change) partial.r(Tal.Or,1:4,cs(gender,age)) #find the partial correlations between the first three variables and the DV (reaction) round(partial.r(Tal.Or[1:4])[4,1:3],2) #The partial correlations with the criterion #Consider the eminence data set from Del Giudice. if(require("psychTools")) { data(eminence) partial.r(reputation ~ works + citations - birth.year, data=eminence) } 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)) phi(c(40,10,10,40)) x <- matrix(c(40,5,20,20),ncol=2) phi(x) 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 #demo <- phi.demo() #compare the phi (lower off diagonal and polychoric correlations # (upper off diagonal) #show the result from tetrachoric which corrects for zero entries by default #round(demo$tetrachoric$rho,2) #show the result from phi2poly #tetrachorics above the diagonal, phi below the diagonal #round(demo$phis,2) cleanEx() nameEx("phi2poly") ### * phi2poly flush(stderr()); flush(stdout()) ### Name: phi2tetra ### Title: Convert a phi coefficient to a tetrachoric correlation ### Aliases: phi2tetra phi2poly ### Keywords: models models ### ** Examples phi2tetra(.3,c(.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 plot.poly plot.irt plot.residuals ### Keywords: multivariate ### ** Examples test.data <- Harman74.cor$cov f4 <- fa(test.data,4) plot(f4) plot(resid(f4)) plot(resid(f4),main="Residuals from a 4 factor solution",qq=FALSE) #not run #data(bfi) #e.irt <- irt.fa(bfi[11:15]) #just the extraversion items #plot(e.irt) #the information curves # ic <- iclust(test.data,3) #shows hierarchical structure plot(ic) #plots loadings # 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 <- fa(circ.data,2) circ.polar <- round(polar(circ.fa),2) circ.polar #compare to the graphic cluster.plot(circ.fa) 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() #compare the phi (lower off diagonal and polychoric correlations (upper off diagonal) #show the result from poly.mat #round(demo$tetrachoric$rho,2) #show the result from phi2poly #tetrachorics above the diagonal, phi below the diagonal #round(demo$phis,2) cleanEx() nameEx("predict.psych") ### * predict.psych flush(stderr()); flush(stdout()) ### Name: predict.psych ### Title: Prediction function for factor analysis, principal components ### (pca), bestScales ### Aliases: predict.psych ### Keywords: multivariate models ### ** Examples set.seed(42) x <- sim.item(12,500) f2 <- fa(x[1:250,],2,scores="regression") # a two factor solution p2 <- principal(x[1:250,],2,scores=TRUE) # a two component solution round(cor(f2$scores,p2$scores),2) #correlate the components and factors from the A set #find the predicted scores (The B set) pf2 <- predict(f2,x[251:500,],x[1:250,]) #use the original data for standardization values pp2 <- predict(p2,x[251:500,],x[1:250,]) #standardized based upon the first set round(cor(pf2,pp2),2) #find the correlations in the B set #test how well these predicted scores match the factor scores from the second set fp2 <- fa(x[251:500,],2,scores=TRUE) round(cor(fp2$scores,pf2),2) pf2.n <- predict(f2,x[251:500,]) #Standardized based upon the new data set round(cor(fp2$scores,pf2.n)) #predict factors of set two from factors of set 1, factor order is arbitrary #note that the signs of the factors in the second set are arbitrary cleanEx() nameEx("predicted.validity") ### * predicted.validity flush(stderr()); flush(stdout()) ### Name: predicted.validity ### Title: Find the predicted validities of a set of scales based on item ### statistics ### Aliases: predicted.validity item.validity validityItem ### Keywords: multivariate models ### ** Examples pred.bfi <- predicted.validity(psychTools::bfi[,1:25], psychTools::bfi[,26:28], psychTools::bfi.keys) pred.bfi cleanEx() nameEx("principal") ### * principal flush(stderr()); flush(stdout()) ### Name: principal ### Title: Principal components analysis (PCA) ### Aliases: principal pca ### Keywords: multivariate models ### ** Examples #Four principal components of the Harman 24 variable problem #compare to a four factor principal axes solution using factor.congruence pc <- principal(Harman74.cor$cov,4,rotate="varimax") mr <- fa(Harman74.cor$cov,4,rotate="varimax") #minres factor analysis pa <- fa(Harman74.cor$cov,4,rotate="varimax",fm="pa") # principal axis factor analysis round(factor.congruence(list(pc,mr,pa)),2) pc2 <- principal(Harman.5,2,rotate="varimax") pc2 round(cor(Harman.5,pc2$scores),2) #compare these correlations to the loadings #now do it for unstandardized scores, and transform obliquely pc2o <- principal(Harman.5,2,rotate="promax",covar=TRUE) pc2o round(cov(Harman.5,pc2o$scores),2) pc2o$Structure #this matches the covariances with the scores biplot(pc2,main="Biplot of the Harman.5 socio-economic variables",labels=paste0(1:12)) #For comparison with SPSS (contributed by Gottfried Helms) pc2v <- principal(iris[1:4],2,rotate="varimax",normalize=FALSE,eps=1e-14) print(pc2v,digits=7) pc2V <- principal(iris[1:4],2,rotate="Varimax",eps=1e-7) p <- print(pc2V,digits=7) round(p$Vaccounted,2) # the amount of variance accounted for is returned as an object of print 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(psychTools::bfi[1:25])) scores <- score.items(keys,psychTools::bfi[1:25]) scores summary(scores) 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.test(50,r) r.test(30,.4,.6) #test the difference between two independent correlations r.test(103,.4,.5,.1) #Steiger case A of dependent correlations r.test(n=103, r12=.4, r13=.5,r23=.1) #for complicated tests, it is probably better to specify correlations by name r.test(n=103,r12=.5,r34=.6,r13=.7,r23=.5,r14=.5,r24=.8) #steiger Case B ##By default, the precision of p values is 2 decimals #Consider three different precisions shown by varying the requested number of digits r12 = 0.693458895410494 r23 = 0.988475791500198 r13 = 0.695966022434845 print(r.test(n = 5105 , r12 = r12 , r23 = r23 , r13 = r13 )) #probability < 0.1 print(r.test(n = 5105 , r12 = r12, r23 = r23 , r13 = r13 ),digits=4) #p < 0.1001 print(r.test(n = 5105 , r12 = r12, r23 = r23 , r13 = r13 ),digits=8) #p< <0.1000759 #an example of how to compare the elements of two matrices R1 <- lowerCor(psychTools::bfi[1:200,1:5]) #find one set of Correlations R2 <- lowerCor(psychTools::bfi[201:400,1:5]) #and now another set sampled #from the same population test <- r.test(n=200, r12 = R1, r34 = R2) round(lowerUpper(R1,R2,diff=TRUE),digits=2) #show the differences between correlations #lowerMat(test$p) #show the p values of the difference between the two matrices adjusted <- p.adjust(test$p[upper.tri(test$p)]) both <- test$p both[upper.tri(both)] <- adjusted round(both,digits=2) #The lower off diagonal are the raw ps, the upper the adjusted ps cleanEx() nameEx("range.correction") ### * range.correction flush(stderr()); flush(stdout()) ### Name: rangeCorrection ### Title: Correct correlations for restriction of range. (Thorndike Case ### 2) ### Aliases: rangeCorrection ### Keywords: multivariate models ### ** Examples rangeCorrection(.33,100.32,48.19) #example from Revelle (in prep) Chapter 4. cleanEx() nameEx("reliability") ### * reliability flush(stderr()); flush(stdout()) ### Name: reliability ### Title: Reports 7 different estimates of scale reliabity including ### alpha, omega, split half ### Aliases: reliability plot.reliability ### Keywords: multivariate models ### ** Examples reliability(psychTools::ability) #an example of finding reliability for all items rel <- reliability(psychTools::ability.keys,psychTools::ability) #use keys to select scales R <- cor(psychTools::ability,use="pairwise") #find the correlations to test rel.R <- reliability(psychTools::ability.keys,R) #this should be the same as rel plot(rel.R) #versus all and subsets all.equal(rel$result.df,rel.R$result.df ) #should be TRUE reliability(psychTools::bfi.keys,psychTools::bfi) #reliability when items are keyed negative ## Not run: ##D ##D #this takes a few seconds but shows nice graphic displays ##D ##D spi.rel <- reliability(psychTools::spi.keys,psychTools::spi,hist=TRUE) #graph them ##D spi.rel #show them ##D ##D #plot them using plot.reliability ##D plot(spi.rel) #draw the density distrbutions ##D ##D plot(spi.rel,split=FALSE) #don't draw the split half density distribution ##D plot(spi.rel,omega=FALSE) # don't add omega values to the diagram ##D #or do this without the densities ##D ##D #plot the first three values in a dot chart ##D error.dots(spi.rel$result.df[,1],sort=FALSE, xlim=c(.3,1),head=16,tail=16, ##D main = expression(paste(omega[h], ~~~~ alpha,~~~~ omega[t]))) ##D #plot the omega_h values ##D error.dots(spi.rel$result.df[,2],sort=FALSE,pch=2,xlim=c(.3,1),head=16,tail=16, ##D main="",labels="",add=TRUE)#add the alpha values ##D error.dots(spi.rel$result.df[,3],sort=FALSE, xlim=c(.3,1),head=16,tail=16, ##D pch=3,labels="", main="",add=TRUE) #and the omega_t values ##D ##D #or, show the smallest and greatest split half, as well as alpha ##D error.dots(spi.rel$result.df[,4],sort=FALSE, xlim=c(.3,1),head=16,tail=16, ##D main = expression(paste(beta, ~~~~ alpha,~~~~ glb))) ##D error.dots(spi.rel$result.df[,5],sort=FALSE,pch=5,xlim=c(.3,1),head=16,tail=16, ##D main="",labels="",add=TRUE)#add the GLB values ##D error.dots(spi.rel$result.df[,2],sort=FALSE,pch=2,xlim=c(.3,1),head=16,tail=16, ##D main="",labels="",add=TRUE)#add the alpha values ##D ##D ## End(Not run) 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) T1 <- rescale(attitude,seq(0,300,50),seq(10,70,10)) #different means and sigmas describe(T1) cleanEx() nameEx("residuals.psych") ### * residuals.psych flush(stderr()); flush(stdout()) ### Name: residuals.psych ### Title: Extract residuals from various psych objects ### Aliases: residuals.psych resid.psych ### Keywords: multivariate models ### ** Examples f3 <- fa(Thurstone,3) residuals(f3) sum(residuals(f3)^2) #include diagonal sum(residuals(f3,diag=FALSE)^2,na.rm=TRUE) #drop diagonal cleanEx() nameEx("reverse.code") ### * reverse.code flush(stderr()); flush(stdout()) ### Name: reverse.code ### Title: Reverse the coding of selected items prior to scale analysis ### Aliases: reverse.code ### Keywords: multivariate ### ** Examples original <- matrix(sample(6,50,replace=TRUE),10,5) keys <- c(1,1,-1,-1,1) #reverse the 3rd and 4th items new <- reverse.code(keys,original,mini=rep(1,5),maxi=rep(6,5)) original[1:3,] new[1:3,] 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) pairs.panels(sat.act) cleanEx() nameEx("scatter.hist") ### * scatter.hist flush(stderr()); flush(stdout()) ### Name: scatterHist ### Title: Draw a scatter plot with associated X and Y histograms, ### densities and correlation ### Aliases: scatter.hist scatterHist ### Keywords: multivariate hplot ### ** Examples data(sat.act) with(sat.act,scatterHist(SATV,SATQ)) scatterHist(SATV ~ SATQ,data=sat.act) #formula input #or for something a bit more splashy scatter.hist(sat.act[5:6],pch=(19+sat.act$gender),col=c("blue","red")[sat.act$gender],grid=TRUE) #better yet scatterHist(SATV ~ SATQ + gender,data=sat.act) #formula input with a grouping variable 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 if(!require(GPArotation)) { message("I am sorry, you must have GPArotation installed to use schmid.")} else { p.jen <- schmid(jen,digits=2) #use the oblimin rotation 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) #deprecated cleanEx() nameEx("score.irt") ### * score.irt flush(stderr()); flush(stdout()) ### Name: scoreIrt ### Title: Find Item Response Theory (IRT) based scores for dichotomous or ### polytomous items ### Aliases: scoreIrt scoreIrt.1pl scoreIrt.2pl score.irt score.irt.poly ### score.irt.2 irt.stats.like make.irt.stats irt.tau irt.se ### Keywords: multivariate models ### ** Examples cleanEx() nameEx("score.items") ### * score.items flush(stderr()); flush(stdout()) ### Name: scoreItems ### Title: Score item composite scales and find Cronbach's alpha, Guttman ### lambda 6 and item whole correlations ### Aliases: scoreItems scoreFast scoreVeryFast score.items ### response.frequencies responseFrequency ### Keywords: multivariate models ### ** Examples #see the example including the bfi data set data(psychTools::bfi) keys.list <- list(agree=c("-A1","A2","A3","A4","A5"), conscientious=c("C1","C2","C3","-C4","-C5"),extraversion=c("-E1","-E2","E3","E4","E5"), neuroticism=c("N1","N2","N3","N4","N5"), openness = c("O1","-O2","O3","O4","-O5")) keys <- make.keys(psychTools::bfi,keys.list) #no longer necessary scores <- scoreItems(keys,psychTools::bfi,min=1,max=6) #using a keys matrix scores <- scoreItems(keys.list,psychTools::bfi,min=1,max=6) # or just use the keys.list summary(scores) #to get the response frequencies, we need to not use the age variable scores <- scoreItems(keys[1:25,],psychTools::bfi[1:25]) #we do not need to specify min or #max if there are no values (such as age) outside the normal item range. scores #The scores themselves are available in the scores$scores object. I.e., describe(scores$scores) #compare this output to that for the impute="none" option for SAPA type data #first make many of the items missing in a missing pattern way missing.bfi <- psychTools::bfi missing.bfi[1:1000,3:8] <- NA missing.bfi[1001:2000,c(1:2,9:10)] <- NA scores <- scoreItems(keys.list,missing.bfi,impute="none",min=1,max=6) scores describe(scores$scores) #the actual scores themselves #If we want to delete scales scores for people who did not answer some items for one #(or more) scales, we can do the following: scores <- scoreItems(keys.list,missing.bfi,totals=TRUE,min=1,max=6) #find total scores describe(scores$scores) #note that missing data were replaced with median for the item scores$scores[scores$missing > 0] <- NA #get rid of cases with missing data describe(scores$scores) 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(psychTools::iqitems) iq.keys <- c(4,4,4, 6,6,3,4,4, 5,2,2,4, 3,2,6,7) score.multiple.choice(iq.keys,psychTools::iqitems) #just convert the items to true or false iq.tf <- score.multiple.choice(iq.keys,psychTools::iqitems,score=FALSE) describe(iq.tf) #compare to previous results cleanEx() nameEx("scoreWtd") ### * scoreWtd flush(stderr()); flush(stdout()) ### Name: scoreWtd ### Title: Score items using regression or correlation based weights ### Aliases: scoreWtd ### Keywords: multivariate models ### ** Examples #find the weights from a regression model and then apply them to a new set #derivation of weights from the first 20 cases model.lm <- lm(rating ~ complaints + privileges + learning,data=attitude[1:20,]) #or use setCor to find the coefficents model <- setCor(rating ~ complaints + privileges +learning,data=attitude[1:20,],std=FALSE) #Apply these to a different set of data (the last 10 cases) #note that the regression coefficients need to be a matrix scores.lm <- scoreWtd(as.matrix(model.lm$coefficients),attitude[21:30,],sums=TRUE,std=FALSE) scores <- scoreWtd(model$coefficients,attitude[21:30,],sums=TRUE,std=FALSE) describe(scores) cleanEx() nameEx("scrub") ### * scrub flush(stderr()); flush(stdout()) ### Name: scrub ### Title: A utility for basic data cleaning and recoding. Changes values ### outside of minimum and maximum limits to NA. ### Aliases: scrub ### Keywords: multivariate ### ** Examples data(attitude) x <- scrub(attitude,isvalue=55) #make all occurrences of 55 NA x1 <- scrub(attitude, where=c(4,5,6), isvalue =c(30,40,50), newvalue = c(930,940,950)) #will do this for the 4th, 5th, and 6th variables x2 <- scrub(attitude, where=c(4,4,4), isvalue =c(30,40,50), newvalue = c(930,940,950)) #will just do it for the 4th column new <- scrub(attitude,1:3,cuts= c(10,40,50,60,100)) #change many values to fewer #get rid of a complicated set of cases and replace with missing values y <- scrub(attitude,where=2:4,min=c(20,30,40),max= c(120,110,100),isvalue= c(32,43,54)) y1 <- scrub(attitude,where="learning",isvalue=55,newvalue=999) #change a column by name y2 <- scrub(attitude,where="learning",min=45,newvalue=999) #change a column by name y3 <- scrub(attitude,where="learning",isvalue=c(45,48), newvalue=999) #change a column by name look for multiple values in that column y4 <- scrub(attitude,where="learning",isvalue=c(45,48), newvalue= c(999,-999)) #change values in one column to one of two different things cleanEx() nameEx("set.cor") ### * set.cor flush(stderr()); flush(stdout()) ### Name: setCor ### Title: Multiple Regression and Set Correlation from matrix or raw input ### Aliases: setCor setCor.diagram set.cor mat.regress matReg ### crossValidation matPlot ### Keywords: models multivariate ### ** Examples #First compare to lm using data input summary(lm(rating ~ complaints + privileges, data = attitude)) setCor(rating ~ complaints + privileges, data = attitude, std=FALSE) #do not standardize z.attitude <- data.frame(scale(attitude)) #standardize the data before doing lm summary(lm(rating ~ complaints + privileges, data = z.attitude)) #regressions on z scores setCor(rating ~ complaints + privileges, data = attitude) #by default we standardize and # the results are the same as the standardized lm R <- cor(attitude) #find the correlations #Do the regression on the correlations #Note that these match the regressions on the standard scores of the data setCor(rating ~ complaints + privileges, data =R, n.obs=30) #now, partial out learning and critical setCor(rating ~ complaints + privileges - learning - critical, data =R, n.obs=30) #compare with the full regression: setCor(rating ~ complaints + privileges + learning + critical, data =R, n.obs=30) #Canonical correlations: #The first Kelley data set from Hotelling kelley1 <- structure(c(1, 0.6328, 0.2412, 0.0586, 0.6328, 1, -0.0553, 0.0655, 0.2412, -0.0553, 1, 0.4248, 0.0586, 0.0655, 0.4248, 1), .Dim = c(4L, 4L), .Dimnames = list(c("reading.speed", "reading.power", "math.speed", "math.power"), c("reading.speed", "reading.power", "math.speed", "math.power"))) lowerMat(kelley1) mod1 <- setCor(y = math.speed + math.power ~ reading.speed + reading.power, data = kelley1, n.obs=140) mod1$cancor #Hotelling reports .3945 and .0688 we get 0.39450592 0.06884787 #the second Kelley data from Hotelling kelley <- structure(list(speed = c(1, 0.4248, 0.042, 0.0215, 0.0573), power = c(0.4248, 1, 0.1487, 0.2489, 0.2843), words = c(0.042, 0.1487, 1, 0.6693, 0.4662), symbols = c(0.0215, 0.2489, 0.6693, 1, 0.6915), meaningless = c(0.0573, 0.2843, 0.4662, 0.6915, 1)), .Names = c("speed", "power", "words", "symbols", "meaningless"), class = "data.frame", row.names = c("speed", "power", "words", "symbols", "meaningless")) lowerMat(kelley) setCor(power + speed ~ words + symbols + meaningless,data=kelley) #formula mode #setCor(y= 1:2,x = 3:5,data = kelley) #order of variables input #Hotelling reports canonical correlations of .3073 and .0583 or squared correlations of # 0.09443329 and 0.00339889 vs. our values of cancor = 0.3076 0.0593 with squared values #of 0.0946 0.0035, setCor(y=c(7:9),x=c(1:6),data=Thurstone,n.obs=213) #easier to just list variable #locations if we have long names #now try partialling out some variables set.cor(y=c(7:9),x=c(1:3),z=c(4:6),data=Thurstone) #compare with the previous #compare complete print out with summary printing sc <- setCor(SATV + SATQ ~ gender + education,data=sat.act) # regression from raw data sc summary(sc) setCor(Pedigrees ~ Sentences + Vocabulary - First.Letters - Four.Letter.Words , data=Thurstone) #showing formula input with two covariates #Do some regressions with real data (rather than correlation matrices) setCor(reaction ~ cond + pmi + import, data = Tal.Or) #partial out importance setCor(reaction ~ cond + pmi - import, data = Tal.Or, main="Partial out importance") #compare with using lm by partialling mod1 <- lm(reaction ~ cond + pmi + import, data = Tal.Or) reaction.import <- lm(reaction~import,data=Tal.Or)$resid cond.import <- lm(cond~import,data=Tal.Or)$resid pmi.import <- lm(pmi~import,data=Tal.Or)$resid mod.partial <- lm(reaction.import ~ cond.import + pmi.import) summary(mod.partial) #lm uses raw scores, so set std = FALSE for setCor print(setCor(y = reaction ~ cond + pmi - import, data = Tal.Or,std = FALSE, main = "Partial out importance"),digits=4) #notice that the dfs of the partial approach using lm are 1 more than the setCor dfs #Show how to find quadratic terms sc <- setCor(reaction ~ cond + pmi + I(import^2), data = Tal.Or) sc #pairs.panels(sc$data) #show the SPLOM of the data #Consider an example of a derivation and cross validation sample set.seed(42) ss <- sample(1:2800,1400) model <- setCor(y=26:28,x=1:25,data=bfi[ss,],plot=FALSE) original.fit <- crossValidation(model,bfi[ss,]) #the derivation set cross.fit <- crossValidation(model,bfi[-ss,]) #the cross validation set summary(original.fit) summary(cross.fit) predicted <- predict(model,bfi[-ss,]) cor2(predicted,bfi[-ss,26:28]) cleanEx() nameEx("sim") ### * sim flush(stderr()); flush(stdout()) ### Name: sim ### Title: Functions to simulate psychological/psychometric data. ### Aliases: sim sim.simplex sim.minor ### Keywords: multivariate datagen ### ** Examples simplex <- sim.simplex() #create the default simplex structure lowerMat(simplex) #the correlation matrix #create a congeneric matrix congeneric <- sim.congeneric() lowerMat(congeneric) R <- sim.hierarchical() lowerMat(R) #now simulate categorical items with the hierarchical factor structure. #Let the items be dichotomous with varying item difficulties. marginals = matrix(c(seq(.1,.9,.1),seq(.9,.1,-.1)),byrow=TRUE,nrow=2) X <- sim.poly.mat(R=R,m=marginals,n=1000) lowerCor(X) #show the raw correlations #lowerMat(tetrachoric(X)$rho) # show the tetrachoric correlations (not run) #generate a structure 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) R <- sim.structure(fx,Phi,fy) cor.plot(R$model) #show it graphically simp <- sim.simplex() #show the simplex structure using cor.plot cor.plot(simp,colors=TRUE,main="A simplex structure") #Show a STARS model simp <- sim.simplex(alpha=.8,lambda=.4) #show the simplex structure using cor.plot cor.plot(simp,colors=TRUE,main="State Trait Auto Regressive Simplex" ) dichot.sim <- sim.irt() #simulate 5 dichotomous items poly.sim <- sim.poly(theta=dichot.sim$theta) #simulate 5 polytomous items that correlate #with the dichotomous items 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) pairs.panels(data.df) #show how the design variables are orthogonal # summary(lm(DV~IV1*IV2*IV3,data=data.df)) summary(aov(DV~IV1*IV2*IV3,data=data.df)) set.seed(42) #demonstrate the effect of not centering the data on the regression data.df <- sim.anova(es1=1,es2=.5,es13=1,center=FALSE) # describe(data.df) # #this one is incorrect, because the IVs are not centered summary(lm(DV~IV1*IV2*IV3,data=data.df)) summary(aov(DV~IV1*IV2*IV3,data=data.df)) #compare with the lm model #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)) summary(aov(DV~IV1*IV2*IV3,data=data.df)) 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 f1=fa(test$observed,scores=TRUE) round(cor(f1$scores,test$latent),2) #factor score estimates are correlated with but not equal to the factor scores set.seed(42) #500 responses to 4 discrete items items <- sim.congeneric(N=500,short=FALSE,low=-2,high=2,categorical=TRUE) d4 <- irt.fa(items$observed) #item response analysis of congeneric measures 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 sim.bonds ### 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 round(jensen,2) set.seed(42) #for reproducible results jensen <- sim.hierarchical(n=10000) #use the same gload and fload values, but produce the data #Compare factor scores using the sl model with those that generated the data lowerCor(jensen$theta) #the correlations of the factors fs <- factor.scores(jensen$observed, jensen$sl) #find factor scores from the data lowerCor(fs$scores) #these are now correlated cor2(fs$scores,jensen$theta) #correlation with the generating factors #compare this to a simulation of the bonds model set.seed(42) R <- sim.bonds() R$R #simulate a non-hierarchical structure 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 #do it again with a hierachical structure gload <- matrix(rep(.7,5) ) five.factor.g <- sim.hierarchical(gload,fload,500,TRUE) #create sample data set #compare these two with omega #not run #om.5 <- omega(five.factor$observed,5) #om.5g <- omega(five.factor.g$observed,5) cleanEx() nameEx("sim.irt") ### * sim.irt flush(stderr()); flush(stdout()) ### Name: sim.irt ### Title: Functions to simulate psychological/psychometric data. ### Aliases: sim.irt sim.rasch sim.npl sim.npn sim.poly sim.poly.npl ### sim.poly.npn sim.poly.ideal sim.poly.ideal.npl sim.poly.ideal.npn ### sim.poly.mat ### Keywords: multivariate datagen ### ** Examples simplex <- sim.simplex() #create the default simplex structure lowerMat(simplex) #the correlation matrix #create a congeneric matrix congeneric <- sim.congeneric() lowerMat(congeneric) R <- sim.hierarchical() lowerMat(R) #now simulate categorical items with the hierarchical factor structure. #Let the items be dichotomous with varying item difficulties. marginals = matrix(c(seq(.1,.9,.1),seq(.9,.1,-.1)),byrow=TRUE,nrow=2) X <- sim.poly.mat(R=R,m=marginals,n=1000) lowerCor(X) #show the raw correlations #lowerMat(tetrachoric(X)$rho) # show the tetrachoric correlations (not run) #generate a structure 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) R <- sim.structure(fx,Phi,fy) cor.plot(R$model) #show it graphically simp <- sim.simplex() #show the simplex structure using cor.plot cor.plot(simp,colors=TRUE,main="A simplex structure") #Show a STARS model simp <- sim.simplex(alpha=.8,lambda=.4) #show the simplex structure using cor.plot cor.plot(simp,colors=TRUE,main="State Trait Auto Regressive Simplex" ) dichot.sim <- sim.irt() #simulate 5 dichotomous items poly.sim <- sim.poly(theta=dichot.sim$theta) #simulate 5 polytomous items that correlate #with the dichotomous items cleanEx() nameEx("sim.item") ### * sim.item flush(stderr()); flush(stdout()) ### Name: sim.item ### Title: Generate simulated data structures for circumplex, spherical, or ### simple structure ### Aliases: sim.spherical item.sim sim.item sim.dichot item.dichot ### sim.circ circ.sim con2cat ### Keywords: multivariate datagen ### ** Examples round(cor(circ.sim(nvar=8,nsub=200)),2) plot(fa(circ.sim(16,500),2)$loadings,main="Circumplex Structure") #circumplex structure # # plot(fa(item.sim(16,500),2)$loadings,main="Simple Structure") #simple structure # cluster.plot(fa(item.dichot(16,low=0,high=1),2)) set.seed(42) data <- mnormt::rmnorm(1000, c(0, 0), matrix(c(1, .5, .5, 1), 2, 2)) #continuous data new <- con2cat(data,c(-1.5,-.5,.5,1.5)) #discreet data polychoric(new) #not run #x12 <- sim.item(12,gloading=.6) #f3 <- fa(x12,3,rotate="none") #f3 #observe the general factor #oblimin(f3$loadings[,2:3]) #show the 2nd and 3 factors. #f3 <- fa(x12,3) #now do it with oblimin rotation #f3 # not what one naively expect. cleanEx() nameEx("sim.multilevel") ### * sim.multilevel flush(stderr()); flush(stdout()) ### Name: sim.multilevel ### Title: Simulate multilevel data with specified within group and between ### group correlations ### Aliases: sim.multilevel sim.multi ### Keywords: multivariate models ### ** Examples #First, show a few results from sim.multi x.df <- sim.multi() #the default is 4 subjects for two variables # over 16 days measured 6 times/day #sb <- statsBy(x.df,group ="id",cors=TRUE) #round(sb$within,2) #show the within subject correlations #get some parameters to simulate data(withinBetween) wb.stats <- statsBy(withinBetween,"Group") rwg <- wb.stats$rwg rbg <- wb.stats$rbg eta <- rep(.5,9) #simulate them. Try this again to see how it changes XY <- sim.multilevel(ncases=100,ngroups=10,rwg=rwg,rbg=rbg,eta=eta) lowerCor(XY$wg) #based upon 89 df lowerCor(XY$bg) #based upon 9 df -- cleanEx() nameEx("sim.omega") ### * sim.omega flush(stderr()); flush(stdout()) ### Name: sim.omega ### Title: Further functions to simulate psychological/psychometric data. ### Aliases: sim.omega sim.general sim.parallel ### Keywords: multivariate datagen ### ** Examples #test <- sim.omega() 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: sim.structure sim.structural sim.correlation simCor ### Keywords: multivariate datagen ### ** Examples #First, create a sem like model with a factor model of x and ys with correlation Phi fx <-matrix(c( .9,.8,.6,rep(0,4),.6,.8,-.7),ncol=2) fy <- matrix(c(.6,.5,.4),ncol=1) rownames(fx) <- c("V","Q","A","nach","Anx") rownames(fy)<- c("gpa","Pre","MA") Phi <-matrix( c(1,0,.7,.0,1,.7,.7,.7,1),ncol=3) #now create this structure gre.gpa <- sim.structural(fx,Phi,fy) print(gre.gpa,2) #correct for attenuation to see structure #the raw correlations are below the diagonal, the adjusted above round(correct.cor(gre.gpa$model,gre.gpa$reliability),2) #These are the population values, # we can also create a correlation matrix sampled from this population GRE.GPA <- sim.structural(fx,Phi,fy,n=250,raw=FALSE) lowerMat(GRE.GPA$r) #or we can show data sampled from such a population GRE.GPA <- sim.structural(fx,Phi,fy,n=250,raw=TRUE) lowerCor(GRE.GPA$observed) congeneric <- sim.structure(f=c(.9,.8,.7,.6)) # a congeneric model congeneric #now take this correlation matrix as a population value and create samples from it example.congeneric <- sim.correlation(congeneric$model,n=200) #create a sample matrix lowerMat(example.congeneric ) #show the correlation matrix #or create another sample and show the data example.congeneric.data <- simCor(congeneric$model,n=200,data=TRUE) describe(example.congeneric.data) lowerCor(example.congeneric.data ) example.skewed <- simCor(congeneric$model,n=200,vars=c(1,2),data=TRUE,skew="log") describe(example.skewed) 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 circ.sim.plot ### Keywords: multivariate datagen ### ** Examples #not run 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: mardia ### Title: Calculate univariate or multivariate (Mardia's test) skew and ### kurtosis for a vector, matrix, or data.frame ### Aliases: mardia skew kurtosi ### Keywords: multivariate models ### ** Examples round(skew(attitude),2) #type 3 (default) round(kurtosi(attitude),2) #type 3 (default) #for the differences between the three types of skew and kurtosis: round(skew(attitude,type=1),2) #type 1 round(skew(attitude,type=2),2) #type 2 mardia(attitude) x <- matrix(rnorm(1000),ncol=10) describe(x) mardia(x) 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 ### Keywords: multivariate ### ** Examples R <- make.hierarchical() round(smc(R),2) cleanEx() nameEx("spider") ### * spider flush(stderr()); flush(stdout()) ### Name: spider ### Title: Make "radar" or "spider" plots. ### Aliases: spider radar ### Keywords: multivariate hplot ### ** Examples op <- par(mfrow=c(3,2)) spider(y=1,x=2:9,data=Thurstone,connect=FALSE) #a radar plot spider(y=1,x=2:9,data=Thurstone) #same plot as a spider plot spider(y=1:3,x=4:9,data=Thurstone,overlay=TRUE) #make a somewhat oversized plot spider(y=26:28,x=1:25,data=cor(psychTools::bfi,use="pairwise"),fill=TRUE,scale=2) par(op) #another example taken from Lippa (2001, page 193) lippa.df <- structure(list(labels = c("Assured - Dominant", "Gregarious\nExtraverted", "Warm\nAgreeable", "Unassuming\nIngeneous", "Unassured - Submissive", "Aloof\nIntroverted", "Cold\nHearted", "Arrogant\nCalculating" ), pos = c(0.8, 0.85, 0.83, 0.8, 0.75, 0.83, 0.85, 0.85), values = c(0.41, -0.29, -0.53, -0.61, -0.38, 0.14, 0.59, 0.6), delta = c(1.1, 1.2, 1.2, 1.1, 1.1, 1.5, 1.2, 1.1)), row.names = c(NA, -8L), class = "data.frame") radar(lippa.df$values,abs=TRUE,labels=lippa.df$labels,angle=90,clockwise=TRUE,lwd=3, label.pos=lippa.df$pos,main="Data from Lippa (2001)",scale=.9,circles=FALSE, cut=0,delta=lippa.df$delta) segments(-1,0,1,0,lwd=.2) # Add hairline axes segments(0,-1,0,1,lwd=.2) text(0,1.05,expression(italic("Masculine Instrumentality"))) text(1.05,0,expression(italic("Feminine Communion")),srt=270) #show how to draw a hexagon RIASEC.df <- structure(list(labels = c("Realistic", "Investigative", "Artistic", "Social", "Enterprising", "Conventional"), Su = c(0.84, 0.26, -0.35, -0.68, 0.04, -0.33), Morris = c(1.14, 0.32, -0.19, -0.38, 0.22, 0.23)), row.names = c(NA, -6L), class = "data.frame") radar(RIASEC.df$Morris,RIASEC.df$labels,clockwise=TRUE,angle=0,absolute=TRUE,circl=FALSE,scale=.7, position=c(1,0,0,0,0,0), lwd=4,label.pos=rep(.80,6),main="",cut=0, shape=TRUE, delta =c(1.1,1.25,1.25, 1.25, 1.45,1.45) ) text(-1.04,0,expression(italic("People")),srt=90) text(1.04,0,expression(italic("Things")),srt=270) text(0,.91,expression(italic("Data"))) text(0,-.91 ,expression(italic("Ideas"))) segments(-1,0,1,0,lwd=.2) #add hairline axes segments(0,-.86,0,.86,lwd=.2) text(0,1.2, "Data from Su") graphics::par(get("par.postscript", pos = 'CheckExEnv')) cleanEx() nameEx("statsBy") ### * statsBy flush(stderr()); flush(stdout()) ### Name: statsBy ### Title: Find statistics (including correlations) within and between ### groups for basic multilevel analyses ### Aliases: statsBy statsBy.boot statsBy.boot.summary faBy ### Keywords: multivariate models ### ** Examples #Taken from Pedhazur, 1997 pedhazur <- structure(list(Group = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), X = c(5L, 2L, 4L, 6L, 3L, 8L, 5L, 7L, 9L, 6L), Y = 1:10), .Names = c("Group", "X", "Y"), class = "data.frame", row.names = c(NA, -10L)) pedhazur ped.stats <- statsBy(pedhazur,"Group") ped.stats #Now do this for the sat.act data set sat.stats <- statsBy(sat.act,c("education","gender"),cors=TRUE) #group by two grouping variables print(sat.stats,short=FALSE) lowerMat(sat.stats$pbg) #get the probability values #show means by groups round(sat.stats$mean) #Do separate factor analyses for each group sb.bfi <- statsBy(psychTools::bfi[1:10], group=psychTools::bfi$gender,cors=TRUE) faBy(sb.bfi,1) #one factor per group faBy(sb.bfi,2) #two factors per group cleanEx() nameEx("structure.diagram") ### * structure.diagram flush(stderr()); flush(stdout()) ### Name: structure.diagram ### Title: Draw a structural equation model specified by two measurement ### models and a structural model ### Aliases: structure.diagram structure.graph structure.sem lavaan.diagram ### sem.diagram sem.graph ### Keywords: multivariate hplot ### ** Examples #A set of measurement and structural models #First set up the various matrices fx <- matrix(c(.9,.8,.7,rep(0,9), .6,.7,-.8,rep(0,9),.5,.6,.4),ncol=3) fy <- matrix(c(.9,.8,.6,rep(0,4),.6,.8,-.7),ncol=2) Phi <- matrix(c(1,.35,0,0,0, .35,1,.5,0,0, 0,.5, 1,0,0, .7,-.6, 0, 1,0, .0, 0, .4,0,1 ),ncol=5,byrow=TRUE) #now draw a number of models f1 <- structure.diagram(fx,main = "A measurement model for x") f2 <- structure.diagram(fx,Phi, main = "A measurement model for x") f3 <- structure.diagram(fy=fy, main = "A measurement model for y") f4 <- structure.diagram(fx,Phi,fy,main="A structural path diagram") f5 <- structure.diagram(fx,Phi,fy,main="A structural path diagram",errors=TRUE) #a mimic model fy <- matrix(c(.9,.8,.6,rep(0,4),.6,.8,-.7),ncol=2) fx <- matrix(c(.6,.5,0,.4),ncol=2) mimic <- structure.diagram(fx,fy=fy,simple=FALSE,errors=TRUE, main="A mimic diagram") fy <- matrix(c(rep(.9,8),rep(0,16),rep(.8,8)),ncol=2) structure.diagram(fx,fy=fy) #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" f2 <- structure.diagram(X2,Phi=phi2,errors=FALSE,main="A symbolic model") #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" f3 <- structure.diagram(X2,Phi=phi2,main="an alternative representation",e.size=.4) #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") f4 <- structure.diagram(X6,phi21,Y3) ###the following example is not run but is included to show how to work with lavaan # and finally, a regression model X7 <- matrix(c("a","b","c","d","e","f"),nrow=6) f5 <- structure.diagram(X7,regression=TRUE,main = "Regression model") #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) f6<- structure.diagram(x8,Phi=r8,regression=TRUE,main="Regression model") 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 phi fy cleanEx() nameEx("super.matrix") ### * super.matrix flush(stderr()); flush(stdout()) ### Name: superMatrix ### Title: Form a super matrix from two sub matrices. ### Aliases: superMatrix superCor 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="") mxy <- superMatrix(mx,my) #show the use of a list to do this as well key1 <- make.keys(6,list(first=c(1,-2,3),second=4:6,all=1:6)) #make a scoring key key2 <- make.keys(4,list(EA=c(1,2),TA=c(3,4))) superMatrix(list(key1,key2)) r <- cor(bfi[1:15],use="pairwise") bfi.scores <- scoreOverlap(bfi.keys[1:2], r,select=FALSE) #note the select = FALSE R <- superCor(bfi.scores,r) lowerMat(R) #or to just get the scale correlations with the items R <- superCor(bfi.scores) round(R,2) 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(psychTools::cubits,labs=c("height","cubit")) describe(cubit) ellipses(cubit,n=1) data(bock) responses <- table2df(bock.table[,2:6],count=bock.table[,7],labs= paste("lsat6.",1:5,sep="")) describe(responses) cleanEx() nameEx("tal_or") ### * tal_or flush(stderr()); flush(stdout()) ### Name: Tal_Or ### Title: Data set testing causal direction in presumed media influence ### Aliases: Tal_Or Tal.Or pmi tctg ### Keywords: datasets ### ** Examples data(Tal.Or) mediate(reaction ~ cond + (pmi), data =Tal.Or,n.iter=50) cleanEx() nameEx("test.irt") ### * test.irt flush(stderr()); flush(stdout()) ### Name: test.irt ### Title: A simple demonstration (and test) of various IRT scoring ### algorthims. ### Aliases: test.irt ### Keywords: multivariate models ### ** Examples #not run #test.irt(9,1000) 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() #not run #test.psych(all=TRUE) # f3 <- fa(bfi[1:15],3,n.iter=5) # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="Varimax") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="varimax") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="bifactor") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="varimin") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="bentlerT") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="geominT") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="equamax") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="Promax") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="cluster") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="biquartimin") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="equamax") # f3 <- fa(bfi[1:15],3,n.iter=5,rotate="Promax") # # fpoly <- fa(bfi[1:10],2,n.iter=5,cor="poly") # f1 <- fa(psychTools::ability,n.iter=4) # f1p <- fa(psychTools::ability,n.iter=4,cor="tet") cleanEx() nameEx("testRetest") ### * testRetest flush(stderr()); flush(stdout()) ### Name: testRetest ### Title: Find various test-retest statistics, including test, person and ### item reliability ### Aliases: testRetest testReliability ### Keywords: multivariate models ### ** Examples cleanEx() nameEx("tetrachor") ### * tetrachor flush(stderr()); flush(stdout()) ### Name: tetrachoric ### Title: Tetrachoric, polychoric, biserial and polyserial correlations ### from various types of input ### Aliases: tetrachoric tetrachor polychoric biserial polydi polyserial ### poly.mat ### Keywords: multivariate ### ** Examples #if(require(mnormt)) { data(bock) tetrachoric(lsat6) polychoric(lsat6) #values should be the same tetrachoric(matrix(c(44268,193,14,0),2,2)) #MPLUS reports.24 #Do not apply continuity correction -- compare with previous analysis! tetrachoric(matrix(c(44268,193,14,0),2,2),correct=0) #the default is to add correct=.5 to 0 cells tetrachoric(matrix(c(61661,1610,85,20),2,2)) #Mplus reports .35 tetrachoric(matrix(c(62503,105,768,0),2,2)) #Mplus reports -.10 tetrachoric(matrix(c(24875,265,47,0),2,2)) #Mplus reports 0 polychoric(matrix(c(61661,1610,85,20),2,2)) #Mplus reports .35 polychoric(matrix(c(62503,105,768,0),2,2)) #Mplus reports -.10 polychoric(matrix(c(24875,265,47,0),2,2)) #Mplus reports 0 #Do not apply continuity correction- compare with previous analysis tetrachoric(matrix(c(24875,265,47,0),2,2), correct=0) polychoric(matrix(c(24875,265,47,0),2,2), correct=0) #the same result #examples from Kirk 1973 #note that Kirk's tables have joint probability followed by marginals, but #tetrachoric needs marginals followed by joint probability tetrachoric(c(.5,.5,.333333)) #should be .5 tetrachoric(c(.5,.5,.1150267)) #should be -.75 tetrachoric(c(.5,.5,.397584)) #should e .8 tetrachoric(c(.158655254,.158655254,.145003)) #should be .99 #the example from Olsson, 1979 x <- as.table(matrix(c(13,69,41,6,113,132,0,22,104),3,3)) polychoric(x,correct=FALSE) #Olsson reports rho = .49, tau row = -1.77, -.14 and tau col = -.69, .67 #give a vector of two marginals and the comorbidity tetrachoric(c(.2, .15, .1)) tetrachoric(c(.2, .1001, .1)) #} else { # message("Sorry, you must have mnormt installed")} # 4 plots comparing biserial to point biserial and latent Pearson correlation set.seed(42) x.4 <- sim.congeneric(loads =c(.9,.6,.3,0),N=1000,short=FALSE) y <- x.4$latent[,1] for(i in 1:4) { x <- x.4$observed[,i] r <- round(cor(x,y),1) ylow <- y[x<= 0] yhigh <- y[x > 0] yc <- c(ylow,yhigh) rpb <- round(cor((x>=0),y),2) rbis <- round(biserial(y,(x>=0)),2) ellipses(x,y,ylim=c(-3,3),xlim=c(-4,3),pch=21 - (x>0), main =paste("r = ",r,"rpb = ",rpb,"rbis =",rbis)) dlow <- density(ylow) dhigh <- density(yhigh) points(dlow$y*5-4,dlow$x,typ="l",lty="dashed") lines(dhigh$y*5-4,dhigh$x,typ="l") } #show non-symmeteric results test1 <- tetrachoric(psychTools::ability[,1:4],psychTools::ability[,5:10]) test2 <- polychoric(psychTools::ability[,1:4],psychTools::ability[,5:10]) all <- tetrachoric(psychTools::ability[,1:10]) cleanEx() nameEx("thurstone") ### * thurstone flush(stderr()); flush(stdout()) ### Name: thurstone ### Title: Thurstone Case V scaling ### Aliases: thurstone ### Keywords: models ### ** Examples data(psychTools::vegetables) thurstone(psychTools::veg) #But consider the case of 100 random orders set.seed((42)) ranks <- matrix(NA,nrow=100,ncol=5) for(i in 1:100) ranks[i,] <- sample(5,5) t <- thurstone(ranks,TRUE) t #show the fits t$hoice #show the choice matrix 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 tr(m) cleanEx() nameEx("unidim") ### * unidim flush(stderr()); flush(stdout()) ### Name: unidim ### Title: Several indices of the unidimensionality of a set of variables. ### Aliases: unidim ### Keywords: models multivariate ### ** Examples #test the unidimensionality of the five factors of the bfi data set. unidim(psychTools::bfi,psychTools::bfi.keys) unidim(psychTools::ability,psychTools::ability.keys) #Try a known 3 factor structure x <- sim.minor(nfact=3,bipolar=FALSE) #this makes all the items positive unidim(x$model) keys.list <- list(first =c(1:4),second = 5:8,third=9:12,all=1:12) unidim(x$model,keys.list) x <- sim.minor(nfact=3) unidim(x$model,keys.list) #we flip the negative items #what about a hierarchical model? H <- sim.hierarchical() # by default, a nice hierarchical model H.keys <- list(First = paste0("V",1:3),Second=paste0("V",4:6),Third=paste0("V",7:9), All = paste0("V",1:9)) unidim(H,H.keys) cleanEx() nameEx("winsor") ### * winsor flush(stderr()); flush(stdout()) ### Name: winsor ### Title: Find the Winsorized scores, means, sds or variances for a ### vector, matrix, or data.frame ### Aliases: winsor winsor.mean winsor.means winsor.sd winsor.var ### Keywords: univar ### ** Examples data(sat.act) winsor.means(sat.act) #compare with the means of the winsorized scores y <- winsor(sat.act) describe(y) xy <- data.frame(sat.act,y) #pairs.panels(xy) #to see the effect of winsorizing x <- matrix(1:100,ncol=5) winsor(x) winsor.means(x) y <- 1:11 winsor(y,trim=.5) cleanEx() nameEx("withinBetween") ### * withinBetween flush(stderr()); flush(stdout()) ### Name: withinBetween ### Title: An example of the distinction between within group and between ### group correlations ### Aliases: withinBetween ### Keywords: datasets ### ** Examples data(withinBetween) pairs.panels(withinBetween,bg=c("red","blue","white","black")[withinBetween[,1]], pch=21,ellipses=FALSE,lm=TRUE) stats <- statsBy(withinBetween,'Group') print(stats,short=FALSE) ### *