https://github.com/cran/nFactors
Tip revision: 592b098fc786911733da1c1953e58c9d1c2e9517 authored by Gilles Raiche on 10 April 2010, 00:00:00 UTC
version 2.3.3
version 2.3.3
Tip revision: 592b098
bentlerParameters.r
bentlerParameters <-
function(x, N, nFactors, log=TRUE, cor=TRUE,
minPar=c(min(lambda) - abs(min(lambda)) +.001, 0.001),
maxPar=c(max(lambda), lm(lambda ~ I(length(lambda):1))$coef[2]),
resParx=c(0.01, 2), resPary=c(0.01, 2),
graphic=TRUE, resolution=30, typePlot="wireframe", ...){
stopMessage <- paste("\n These indices are only valid with a principal component solution.\n",
" ...................... So, only positive eugenvalues are permitted.\n",
sep="")
lambda <- eigenComputes(x, cor=cor, ...)
if (length(which(lambda <0 )) > 0) {cat(stopMessage);stop()}
k <- nFactors
p <- length(lambda)
q <- p-k
i <- 1:q
x <- q-i
l <- lambda[k+i]
n <- N - 1
# Bentler (1996, p. 133) maximization of equations 8 and 9
f1 <- function(n,l,x,alpha,beta) sum((n*l-(n+1)*(alpha+beta*x))/((alpha+beta*x)^2))
f2 <- function(n,l,x,alpha,beta) sum((n*l-(n+1)*(alpha+beta*x))*x/((alpha+beta*x)^2))
f <- function(alpha,beta) f1(n,l,x,alpha,beta)^2+f2(n,l,x,alpha,beta)^2
if (log == FALSE) F <- function(y) f(y[1],y[2]) else F <- function(y) log(f(y[1],y[2]))
figure <- NULL
if (graphic == TRUE) {
p1 <- seq(resParx[1], resParx[2], length=resolution)
p2 <- seq(resPary[1], resPary[2], length=resolution)
data <- expand.grid(Alpha = p1, Beta = p2)
data <- data.frame(data, y=numeric(length(data$Alpha)))
for( i in 1:length(data$Alpha)) data$y[i] <- F(c(data$Alpha[i],data$Beta[i]))
if (log == FALSE) zlab <- "y" else zlab <- "log(y)"
if (typePlot == "wireframe") figure <- wireframe( y ~ Alpha * Beta, data=data, zlab=zlab, ...)
if (typePlot == "contourplot") figure <- contourplot(y ~ Alpha * Beta, data=data, region=TRUE, ...)
if (typePlot == "levelplot") figure <- levelplot( y ~ Alpha * Beta, data=data, region=TRUE, ...)
}
res <- nlminb(objective=F,start=lm(l~x)$coefficients,lower=c(minPar[1],minPar[2]),upper=c(maxPar[1],maxPar[2]))
para <- res$par[1]
parb <- res$par[2]
# Bentler (1996, p. 133) equation 7
# !!! Warning: Bentler and Yuan (1998) were in error for the definition of LRT !!!
# !!! So N and n must be inversed in the first logarithm !!!
lrt <- N*(k-p)*(log(n/N)+1)-N*sum(log(lambda[(k+1):p]/(para+parb*x))) + n*sum(lambda[(k+1):p]/(para+parb*x))
df <- q-2
resp <- list(convergence=res$convergence, figure=figure, coefficients=res$par,
lrt=lrt, df=df,k=k,p.value=1-pchisq(lrt,df))
names(resp$coefficients)<-c("alpha","beta")
return(resp)
}