R/fCFA.R In cfa: Configural Frequency Analysis (CFA)

Defines functions `fCFA`

````fCFA` <-
function(m.i, X, tabdim, alpha = 0.05)
{

#tabdim... a vector with table dimensions, e.g., tabdim = c(3,2,3)

#if (alpha.cor == "Bonferroni") alpha <- alpha/length(m.i)

resList <- NULL
chisqVec <- NULL
chidfVec <- NULL
pvalueVec <- NULL
strucMat <- NULL
devVec <- NULL
typevec <- NULL
i <- -1

fit <- FALSE
startdim <- dim(X)[2]
res <- "Pearson"
#bonf <- length(m.i)
while (fit==FALSE)
{
i <- i+1
result <- glm.fit(X,m.i,family=poisson())              #model fit
efreq <- result\$fitted.values                            #expected frequencies

if (res == "Pearson") {
stres <- (m.i-efreq)/sqrt(efreq)                     #standardized Pearson residuals
stresa <- abs(stres)                                   #standardized residuals (absolute)
pvec <- 1-pnorm(stresa)
}
##if (res == "Poisson") {                                  #standardized Poisson residuals
##      w.i <- result\$weights
##      W <- diag(w.i)
##      Hat <- (mroot(W))%*%X%*%(solve(t(X)%*%W%*%X))%*%t(X)%*%(mroot(W))
##      h.i <- diag(Hat)
##      stres <- (m.i-efreq)/sqrt(efreq*(1-h.i))
##      stresa <- abs(stres)
##      pvec <- 1-pnorm(stresa)
##    }
##    if (res == "Deviance") {
##      d.i <- 2*m.i*log(m.i/efreq)
##      stres <- sqrt(abs(d.i))*sign(m.i-efreq)
##      stresa <- abs(stres)
##      pvec <- 1-pnorm(stresa)
##    }

chisq <- sum((m.i-efreq)^2/efreq)                      #chi-square value
chisqVec <- cbind(chisqVec,chisq)
devVec <- cbind(devVec,result\$deviance)
chidfVec <- cbind(chidfVec,result\$df.residual)
pvalue <- 1-pchisq(result\$deviance,result\$df.residual)                          #pvalue
pvalueVec <- cbind(pvalueVec,pvalue)
fitvec <- c(chisq,result\$deviance,result\$df.residual,pvalue)
names(fitvec) <- c("LR","X^2","df","p")
Xdes <- as.matrix(X)
colnames(Xdes) <- NULL
rownames(Xdes) <- NULL
rL <- list(desmat = Xdes, exp.freq = efreq, fitvec = fitvec)
resList <- c(resList,rL)

if (pvalue < alpha) {
strucvec <- as.integer(stresa==max(stresa))            #cells to be blanked out
strucMat <- rbind(strucMat,strucvec)
X <- cbind(X,strucvec)                                 #new design matrix
fit <- FALSE
if (m.i[strucvec==1] > efreq[strucvec==1]) {
typ = "type"
} else {
typ = "antitype"
}
typevec <- c(typevec,typ)
} else {
fit <- TRUE
}
if (result\$df.residual==0) {
fit <- TRUE                   #no more df left
warning("No more df left. Iteration is abandoned.")
}
}

if (is.null(strucMat)) stop("Base model fits! No types/antitypes found!\n",call.=FALSE)

rownames(strucMat) <- paste("Step",1:(dim(strucMat)[1]))     #matrix with excluded elements
td <- length(tabdim)
gridlist <- tapply(tabdim[td:1],1:td,function(x) 1:x)
namesmat <- expand.grid(gridlist)[,td:1]
strvec <- apply(namesmat,1,paste,collapse="")
colnames(strucMat) <- strvec

restable <- data.frame(as.vector(devVec), as.vector(chisqVec), as.vector(chidfVec), as.vector(pvalueVec))
dimnames(restable)[[2]] <- c("LR","X^2","df","p")
dimnames(restable)[[1]] <- paste("Step",0:(dim(restable)[1]-1))
names(typevec) <- paste("Step",1:length(typevec))

desmat <- resList[(length(resList)-2):length(resList)][[1]]               #design matrix

lcomp <- unique(names(resList))
nsteps <- dim(restable)[1]
steplab <- paste("step", 1:nsteps, sep = "")
names(resList) <- as.vector(t(outer(steplab, lcomp, paste, sep = "")))

result <- list(restable = restable, design.mat = desmat,  struc.mat = strucMat,
typevec = typevec, resstep = resList)

class(result) <- "fCFA"

result
}
```

Try the cfa package in your browser

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

cfa documentation built on May 2, 2019, 1:46 p.m.