/
omega.R
201 lines (175 loc) · 11 KB
/
omega.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
"omegah" <-
function(m,nfactors=3,fm="minres",key=NULL,flip=TRUE, digits=2,title="Omega",sl=TRUE,labels=NULL, plot=TRUE,n.obs=NA,rotate="oblimin",Phi = NULL,option="equal",covar=FALSE,two.ok=FALSE,...) {
#m is a correlation matrix, or if not, the correlation matrix is found
#nfactors is the number of factors to extract
#key allows items to be reversed scored if desired
#if Phi is not null, this implies that we have been given a factor matrix -- added May 30, 2010
if(!requireNamespace('GPArotation') && (rotate !="cluster")) {stop("I am sorry, you need to have the GPArotation package installed")}
cl <- match.call()
nvar <- dim(m)[2]
raw.data <- NULL
if(is.null(Phi)) { #the normal case is to do the factor analysis of the raw data or the correlation matrix
if(!isCorrelation(m)) { #should also check for covariance matrix
n.obs <- dim(m)[1]
m <- as.matrix(m)
raw.data <- m #added 9/1/14
if(covar) {m <- cov(m,use="pairwise") } else {m <- cor(m,use="pairwise")}
} else {
if(!covar) m <- cov2cor(as.matrix(m)) #make sure it is a correlation matrix not a covariance or data matrix (if we change this, we will need to change the calculation for omega later)
}
if(is.null(colnames(m))) { rownames(m) <- colnames(m) <- paste("V",1:nvar,sep="") }
m.names <- colnames(m)
if (!is.null(key)) { m <- diag(key) %*% m %*% diag(key)
colnames(m) <- m.names #flip items if we choose to do so
flip <- FALSE #we do this if we specify the key
} else {key <- rep(1,nvar) }
signkey <- strtrim(key,1)
signkey[signkey=="1"] <- ""
m.names <- paste(m.names,signkey,sep="")
colnames(m) <- rownames(m) <- m.names
if ((nvar < 6) && (fm =="mle") ) {message(paste("In omega, 3 factors are too many for ",nvar," variables using mle. Using minres instead",sep=""))
fm <- "minres"}
} else { m.names <- rownames(m) } #add the names if we have a factor input
gf <-schmid(m,nfactors,fm,digits,rotate=rotate,n.obs=n.obs,Phi=Phi,option=option,covar=covar,two.ok=two.ok, ...)
if(!is.null(Phi)) { model <- m
nfactors <- dim(model)[2]
m <- factor.model(model,Phi=Phi,U2=FALSE) #estimate the correlation matrix from the factor model
nvar <- dim(m)[2]
if(is.null(rownames(m))) {colnames(m) <- rownames(m) <- paste("V",1:nvar)}
}
gload <- gf$sl[,1]
if (flip) { #should we think about flipping items ?
key <- sign(gload)
key[key==0] <- 1 # a rare and weird case where the gloading is 0 and thus needs not be flipped
if (sum(key) < nvar) { #some items have negative g loadings and should be flipped
m <- diag(key) %*% m %*% diag(key) #this is just flipping the correlation matrix so we can calculate alpha
gf$sl[,1:(nfactors+1)] <- diag(key) %*% gf$sl[,1:(nfactors+1)]
signkey <- strtrim(key,1)
signkey[signkey=="1"] <- ""
m.names <- paste(m.names,signkey,sep="")
colnames(m) <- rownames(m) <- m.names
rownames(gf$sl) <- m.names
}
}
Vt <- sum(m) #find the total variance in the scale
Vitem <- sum(diag(m))
gload <- gf$sl[,1]
gsq <- (sum(gload))^2
uniq <- sum(gf$sl[,"u2"])
if((nfactors == 1) && (fm=="pc")) {gsq <- Vt - uniq
warning("omega_h is not meaningful for a principal components analysis with one component")} #weird condition when using fm=pc and 1 factor
om.tot <- (Vt-uniq)/Vt
om.limit <- gsq/(Vt-uniq)
alpha <- ((Vt-Vitem)/Vt)*(nvar/(nvar-1))
sum.smc <- sum(smc(m,covar=covar))
lambda.6 <- (Vt +sum.smc-sum(diag(m)))/Vt
if (!is.null(digits)) {omega <-list(omega_h= gsq/Vt,alpha=alpha,lambda.6 = lambda.6,omega.tot =om.tot,schmid=gf ,key = key,title=title)
dg <-max(digits-1,1)} else {
omega <- list(omega_h= gsq/Vt,alpha=alpha,omega.tot=om.tot,schmid=gf,key=key,title=title)
dg <- 1}
ev <- colSums(gf$sl[,1:(nfactors+1)]^2)
ECV <- ev[1]/sum(ev)
omega.stats <- factor.stats(m,gf$sl[,1:(nfactors+1)],n.obs=n.obs)
general.stats <- factor.stats(m,as.matrix(gf$sl[,1]),n.obs=n.obs) #just get fit for the general factor
if (nfactors<2) plot <- FALSE
# if(require(Rgraphviz) && plot) {omega.model <-omega.graph(omega,title=title,sl=sl,labels=labels,digits=dg) } else {omega.model <- omega.sem(omega,sl=sl)}
omega.model <- omega.sem(omega,sl=sl)
#find the subset omegas
omg <- omgo <- omt<- rep(NA,nfactors+1)
sub <- apply(gf$sl,1,function(x) which.max(abs(x[2:(nfactors+1)])))
grs <- 0
for(group in( 1:nfactors)) {
groupi <- which(sub==group)
if(length(groupi) > 0) {
Vgr <- sum(m[groupi,groupi])
gr <- sum(gf$sl[groupi,(group+1)])
grs <- grs + gr^2
omg[group+1] <- gr^2/Vgr
omgo[group+1] <- sum(gf$sl[groupi,1])^2/Vgr
omt[group+1] <- (gr^2+ sum(gf$sl[groupi,1])^2)/Vgr
}
omgo[1] <- sum(gf$sl[,1])^2/sum(m) #omega h
omg[1] <- grs/sum(m) #omega of subscales
omt[1] <- om.tot
om.group <- data.frame(total=omt,general=omgo,group=omg)
rownames(om.group) <- colnames(gf$sl)[1:(nfactors+1)]
} #moved after tge bext line (6/21/18)
#we should standardize the raw.data before doing the next step
#change this to call factor.scores
if(!is.null(raw.data)) { #added mean imputation 4/22/23
miss <- which(is.na(raw.data),arr.ind=TRUE)
item.means <- colMeans(raw.data,na.rm=TRUE) #replace missing values with means
raw.data[miss]<- item.means[miss[,2]]
z.data <- scale(raw.data)
scores <- z.data %*% omega.stats$weights
} else {scores<- NULL}
# if(!is.null(raw.data)) {scores <- raw.data %*% omega.stats$weights} else {scores<- NULL}
# }
omega <- list(omega_h= gsq/Vt,omega.lim = om.limit,alpha=alpha,omega.tot=om.tot,G6=lambda.6,schmid=gf,key=key,stats = omega.stats,ECV=ECV,gstats = general.stats,call=cl,title=title,R = m,model=omega.model,omega.group=om.group,scores=scores,Call=cl)
class(omega) <- c("psych","omega")
if(plot) omega.diagram(omega,main=title,sl=sl,labels=labels,digits=dg)
return(omega)
}
#April 4, 2011 added a check for fm=pc and nfactors == 1 to solve problem of omega_h < omega_t -- probably not a good idea. removed
#January 9, 2014 added omega scores if the raw data are given
"omega" <-
function(m,nfactors=3,fm="minres",n.iter=1,p=.05,poly=FALSE,key=NULL,flip=TRUE, digits=2,title="Omega",sl=TRUE,labels=NULL, plot=TRUE,n.obs=NA,rotate="oblimin",Phi = NULL,option="equal",covar=FALSE,...) {
cl <- match.call()
if(is.data.frame(m) || is.matrix(m)) {if((isCorrelation(m)) | (isCovariance(m) && covar)) {if(is.na(n.obs) && (n.iter > 1)) stop("You must specify the number of subjects if giving a correlation matrix")
# if(!require(MASS)) stop("You must have MASS installed to simulate data from a correlation matrix")
}
}
if(!is.data.frame(m) && !is.matrix(m)) { n.obs=m$n.obs
if(poly) {
pol <- list(rho=m$rho,tau = m$tau,n.obs=m$n.obs)
m <- m$rho
} else { m <- m$R
}
} else { #new data
if(poly) { pol <- polychoric(m)
m <- pol$rho
n.obs <- pol$n.obs} }
om <- omegah(m=m,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=n.obs,rotate=rotate,Phi = Phi,option=option,covar=covar,...) #call omega with the appropriate parameters
if(is.na(n.obs) ) {n.obs <- om$stats$n.obs}
replicates <- list()
if(n.iter > 1) {for (trials in 1:n.iter) {
if(dim(m)[1] == dim(m)[2]) {#create data sampled from multivariate normal with correlation
nvar <- dim(m)[1]
#mu <- rep(0, nvar)
# m <- mvrnorm(n = n.obs, mu, Sigma = m, tol = 1e-06, empirical = FALSE)
#the next 3 lines replaces mvrnorm (taken from mvrnorm, but without the checks)
eX <- eigen(m)
m <- matrix(rnorm(nvar * n.obs),n.obs)
m <- t(eX$vectors %*% diag(sqrt(pmax(eX$values, 0)), nvar) %*% t(m))
} else {m <- m[sample(n.obs,n.obs,replace=TRUE),]}
if(poly) {pol <- polychoric(m)
oms <- omegah(m=pol$rho,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=pol$n.obs,rotate=rotate,Phi = Phi,option=option,...) #call omega with the appropriate parameters
} else {
oms <- omegah(m=m,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=n.obs,rotate=rotate,Phi = Phi,option=option,...) #call omega with the appropriate parameters
}
# oms <-omegah(m=m,nfactors=nfactors,fm=fm,key=key,flip=flip, digits=digits,title=title,sl=sl,labels=labels, plot=plot,n.obs=n.obs,rotate=rotate,Phi = Phi,option=option,...) #call fa with the appropriate parameters
replicates[[trials]] <- list(omega=oms$omega_h,alpha=oms$alpha,omega.tot=oms$omega.tot,G6=oms$G6,omega.lim=oms$omega.lim)
}
replicates <- matrix(unlist(replicates),ncol=5,byrow=TRUE)
z.replicates <- cbind(fisherz(replicates[,1:4]),replicates[,5]) #convert to z scores
means <- colMeans(z.replicates,na.rm=TRUE)
sds <- apply(z.replicates,2,sd,na.rm=TRUE)
ci.lower <- means + qnorm(p/2) * sds
ci.upper <- means + qnorm(1-p/2) * sds
ci <- data.frame(lower = ci.lower,upper=ci.upper)
ci <- rbind(fisherz2r(ci[1:4,]),ci[5,])
rownames(ci) <- c("omega_h","alpha","omega_tot","G6","omega_lim")
colnames(replicates) <- names(means) <- names(sds) <- rownames(ci)
conf <- list(means = means,sds = sds,ci = ci,Call= cl,replicates=replicates)
om$Call=cl
results <- list(om = om,ci=conf) } else {om$Call=cl
if(poly) {om$rho <- pol$rho
om$tau <- pol$tau
om$n.obs <- pol$n.obs
}
results <- om}
class(results) <- c("psych","omega")
return(results)
}
#written April 25, 2011
#adapted May 12, 2011 to be the primary version of omega