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
|
# copied from OpenMx R/MxSummary.R
pChiSqFun <- function(x, val, degf, goal){
goal - pchisq(val, degf, ncp=x)
}
rmseaConfidenceIntervalHelper <- function(chi.squared, df, N, lower, upper){
# Lower confidence interval
if( pchisq(chi.squared, df=df, ncp=0) >= upper){ #sic
lower.lam <- uniroot(f=pChiSqFun, interval=c(1e-10, 1e4), val=chi.squared, degf=df, goal=upper)$root
# solve pchisq(ch, df=df, ncp=x) == upper for x
} else{
lower.lam <- 0
}
# Upper confidence interval
if( pchisq(chi.squared, df=df, ncp=0) >= lower){ #sic
upper.lam <- uniroot(f=pChiSqFun, interval=c(1e-10, 1e4), val=chi.squared, degf=df, goal=lower)$root
# solve pchisq(ch, df=df, ncp=x) == lower for x
} else{
upper.lam <- 0
}
lower.rmsea <- sqrt(lower.lam/(N*df))
upper.rmsea <- sqrt(upper.lam/(N*df))
return(c(lower=lower.rmsea, upper=upper.rmsea))
}
computeFitStatistics <- function(likelihood, DoF, numObs,
independence, indDoF, saturated=0, satDoF=0) {
chi <- likelihood - saturated
chiDoF <- DoF - satDoF # DoF = obsStat-model.ep; satDoF = obsStat-sat.ep; So sat.ep-model.ep == DoF-satDoF
CFI <- (independence - indDoF - likelihood + DoF)/(independence - indDoF - saturated + satDoF)
TLI <- 1
rmseaSquared <- 0
RMSEA <- 0
RMSEACI <- c(lower=NA, upper=NA)
if (!is.na(chiDoF) && chiDoF > 0) {
TLI <- ((independence-saturated)/(indDoF-satDoF) - (chi)/(DoF-satDoF))/((independence-saturated)/(indDoF-satDoF) - 1)
# Here we use N in the denominator as given in the original
# RMSEA paper. The difference between N and N-1 is negligible
# for sample sizes over 30. RMSEA should not be taken seriously
# such small samples anyway.
rmseaSquared <- (chi / (chiDoF) - 1) / numObs
RMSEACI <- c(lower=NA, upper=NA)
if (length(rmseaSquared) == 0 || is.na(rmseaSquared) ||
is.nan(rmseaSquared)) {
# || (rmseaSquared < 0)) { # changed so 'rmseaSquared < 0' yields zero with comment
RMSEA <- NA
} else if (rmseaSquared < 0) {
RMSEA <- 0.0
} else {
RMSEA <- sqrt(rmseaSquared)
ci <- try(rmseaConfidenceIntervalHelper(chi, chiDoF, numObs, .025, .975))
if (!inherits(ci, "try-error")) RMSEACI <- ci
}
}
list(CFI=CFI, TLI=TLI, RMSEA=RMSEA, RMSEASquared=rmseaSquared, RMSEACI=RMSEACI)
}
catFitStatistics <- function(x) {
cat("CFI:", x$CFI, '\n')
cat("TLI:", x$TLI, '\n')
if (length(x$RMSEASquared) == 1 && !is.na(x$RMSEASquared) && x$RMSEASquared < 0.0) {
cat("RMSEA: ", x$RMSEA, '*(Non-centrality parameter is negative)', '\n')
} else {
cat("RMSEA: ", x$RMSEA, " [95% CI (", x$RMSEACI[1], ", ", x$RMSEACI[2], ")]", '\n', sep="")
}
}
|
