R/gof.R

In kkholst/lava: Latent Variable Models

##' @export
rsq <- function(x,stderr=FALSE) {

    if (stderr) {
        v <- endogenous(x)
        vpar <- paste(v,v,sep=lava.options()$symbol[2])
        ic.v <- IC(model.frame(x)[,v])
        ic.mod <- IC(x)
        coef0 <- c(attributes(ic.v)$coef[vpar],
                   coef(x)[vpar])
        ic0 <- cbind(ic.v[,vpar],ic.mod[,vpar])
        p <- length(v)
        idx <- seq_len(p);
        ee <- estimate(NULL,data=NULL,
                       function(x) {
                           res <- (x[idx]-x[idx+p])/x[idx]
                           names(res) <- v
                           as.list(res)
                       },
                       print=function(x,...) {
                           cat("\nR-squared:\n\n")
                           print(x$coefmat)
                       },
                       coef=coef0, IC=ic0)

        res <- ee
        ##res <- list(ee)
        ## for (lat in latent(x)) {

        ##     v <- intersect(children(x,lat),endogenous(x))
        ##     vpar <- paste(v,v,sep=lava.options()$symbol[2])
        ##     lpar <- paste(lat,lat,sep=lava.options()$symbol[2])
        ##     rpar <- paste(v,lat,sep=lava.options()$symbol[1])
        ##     fix <- c(x$model$fix[lat,v,drop=TRUE],x$model$covfix[lat,lat])
        ##     pp <- coef(x)
        ##     idx <- x$model$parpos$A[lat,v]
        ##     idx2 <- x$model$parpos$P[lat,lat]
        ##     p0 <- c(idx,idx2)
        ##     p1 <- setdiff(unique(p0),0)
        ##     p2 <- match(p0,p1)

        ##     k <- length(v)
        ##     coef0 <- c(pp[p1],attributes(ic.v)$coef[vpar])
        ##     ic0 <- cbind(ic.mod[,p1],ic.v[,vpar])
        ##     ee <- estimate(NULL,data=NULL,
        ##                    function(p) {
        ##                        p. <- p[p2]
        ##                        p.[is.na(p.)] <- fix[is.na(p.)]
        ##                        res <- p.[seq_len(k)]^2*p.[k+1]/tail(p,k)
        ##                        names(res) <- v
        ##                        as.list(res)
        ##                    },
        ##                    print=function(x,...) {
        ##                        cat("\nVariance explained by '", lat,"':\n\n",sep="")
        ##                        print(x$coefmat)
        ##                    },coef=coef0,IC=ic0)
        ##     res <- c(res,list(ee))
        ## }

        return(res)
    }



    v <- c(endogenous(x),setdiff(latent(x),parameter(Model(x))))
    res <- coef(x,9,std="yx")
    idx <- with(attributes(res),
                which(type=="variance" & (var==from)))
    nam <- attributes(res)$var[idx]
    res <- 1-res[idx,5]
    names(res) <- nam
    res <- list("R-squared"=res)
    ## M <- moments(x,coef(x))
    ## v <- setdiff(vars(x),exogenous(x))
    ## vvar <- M$Cfull[cbind(v,v)]
    ## rsq <- (vvar-M$P[cbind(v,v)])/vvar

    if (length(latent(x))>0) {
        M <- moments(x,coef(x))
        nn <- names(res)
        for (lat in latent(x)) {
            v <- intersect(children(x,lat),endogenous(x))
            varl <- M$Cfull[lat,lat]
            varv <- M$Cfull[cbind(v,v)]
            fix <- c(x$model$fix[lat,v,drop=TRUE])
            pp <- coef(x)
            if (inherits(x,"lvm.missing")) {
                mp <- match(coef(x$model),names(coef(x)))
                pp <- pp[mp]
            }
            idx1 <- x$model$parpos$A[lat,v]
            ##idx2 <- x$model$parpos$P[lat,lat]
            ##idx3 <- x$model$parpos$P[cbind(v,v)]
            p0 <- c(idx1)
            p1 <- setdiff(unique(p0),0)
            p2 <- match(p0,p1)
            p <- pp[p1]
            p. <- p[p2]
            p.[is.na(p.)] <- fix[is.na(p.)]
            k <- length(v)
            val <- (p.^2*varl)/varv; names(val) <- v
            res <- c(res,list(val))
            nn <- c(nn,paste0("Variance explained by '",lat,"'"))
        }
        names(res) <- nn
    }
    res
}

satmodel <- function(object,logLik=TRUE,data=model.frame(object),
              control=list(trace=1),
              weights=Weights(object),estimator=object$estimator,
              missing=inherits(object,"lvm.missing"),
              regr=FALSE,
              ...) {
    if (object$estimator=="gaussian" & logLik & !missing) {
        if (class(object)[1]%in%c("multigroupfit","multigroup")) {

            ll <- structure(0,nall=0,nobs=0,df=0,class="logLik")
            for (i in seq_len(Model(object)$ngroup)) {
                l0 <- logLik(Model(Model(object))[[i]],data=model.frame(object)[[i]],type="sat")

                ll <- ll+l0
                for (atr in c("nall","nobs","df"))
                    attributes(ll)[[atr]] <- attributes(ll)[[atr]]+attributes(l0)[[atr]]
            }

        }
        return(logLik(object, type="sat"))
    }
    covar <- exogenous(object)
    y <- endogenous(object)
    m0 <- Model(object)
    if (length(covar)>0)
        suppressWarnings(m0 <- regression(m0,y,covar))
    if (length(latent(m0))>0)
        kill(m0) <- latent(m0)
    cancel(m0) <- y
    if (!regr)
        suppressWarnings(covariance(m0) <- y)
    else {
        if (length(y)>1) {
            for (i in seq_len(length(y)-1))
                for (j in seq(i+1,length(y))) {
                    m0 <- regression(m0,y[i],y[j])
                }
        }
        exogenous(m0) <- covar
    }
    if (is.null(control$start)) {
        mystart <- rep(0,with(index(m0), npar.mean+npar))
        mystart[variances(m0,mean=TRUE)] <- 1
        control$start <- mystart
    }
    message("Calculating MLE of saturated model:\n")
    e0 <- estimate(m0,data=data,weights=weights,estimator=estimator,messages=0,control=control,missing=missing,...)
    if (logLik)
        return(logLik(e0))
    return(e0)
}

condition <- function(x) {
    a <- svd(x)
    return(max(a$d)/min(a$d))
}


##' Extract model summaries and GOF statistics for model object
##'
##' Calculates various GOF statistics for model object including global
##' chi-squared test statistic and AIC. Extract model-specific mean and variance
##' structure, residuals and various predicitions.
##'
##'
##' @aliases gof gof.lvmfit moments moments.lvm information information.lvmfit
##' score score.lvmfit logLik.lvmfit
##' @param object Model object
##' @param x Model object
##' @param p Parameter vector used to calculate statistics
##' @param data Data.frame to use
##' @param latent If TRUE predictions of latent variables are included in output
##' @param data2 Optional second data.frame (only for censored observations)
##' @param weights Optional weight matrix
##' @param n Number of observations
##' @param conditional If TRUE the conditional moments given the covariates are
##' calculated. Otherwise the joint moments are calculated
##' @param model String defining estimator, e.g. "gaussian" (see
##' \code{estimate})
##' @param debug Debugging only
##' @param chisq Boolean indicating whether to calculate chi-squared
##' goodness-of-fit (always TRUE for estimator='gaussian')
##' @param level Level of confidence limits for RMSEA
##' @param rmsea.threshold Which probability to calculate, Pr(RMSEA<rmsea.treshold)
##' @param all Calculate all (ad hoc) FIT indices: TLI, CFI, NFI, SRMR, ...
##' @param \dots Additional arguments to be passed to the low level functions
##' @usage
##'
##' gof(object, ...)
##'
##' \method{gof}{lvmfit}(object, chisq=FALSE, level=0.90, rmsea.threshold=0.05,all=FALSE,...)
##'
##' moments(x,...)
##'
##' \method{moments}{lvm}(x, p, debug=FALSE, conditional=FALSE, data=NULL, latent=FALSE, ...)
##'
##' \method{logLik}{lvmfit}(object, p=coef(object),
##'                       data=model.frame(object),
##'                       model=object$estimator,
##'                       weights=Weights(object),
##'                       data2=object$data$data2,
##'                           ...)
##'
##' \method{score}{lvmfit}(x, data=model.frame(x), p=pars(x), model=x$estimator,
##'                    weights=Weights(x), data2=x$data$data2, ...)
##'
##' \method{information}{lvmfit}(x,p=pars(x),n=x$data$n,data=model.frame(x),
##'                    model=x$estimator,weights=Weights(x), data2=x$data$data2, ...)
##'
##' @return A \code{htest}-object.
##' @author Klaus K. Holst
##' @keywords methods models
##' @export
##' @examples
##' m <- lvm(list(y~v1+v2+v3+v4,c(v1,v2,v3,v4)~x))
##' set.seed(1)
##' dd <- sim(m,1000)
##' e <- estimate(m, dd)
##' gof(e,all=TRUE,rmsea.threshold=0.05,level=0.9)
##'
##'
##' set.seed(1)
##' m <- lvm(list(c(y1,y2,y3)~u,y1~x)); latent(m) <- ~u
##' regression(m,c(y2,y3)~u) <- "b"
##' d <- sim(m,1000)
##' e <- estimate(m,d)
##' rsq(e)
##' ##'
##' rr <- rsq(e,TRUE)
##' rr
##' estimate(rr,contrast=rbind(c(1,-1,0),c(1,0,-1),c(0,1,-1)))
##'
`gof` <-
    function(object,...) UseMethod("gof")

##' @export
gof.lvmfit <- function(object,chisq=FALSE,level=0.90,rmsea.threshold=0.05,all=FALSE,...) {
    n <- object$data$n
    if (class(object)[1]=="multigroupfit") n <- sum(unlist(lapply(object$model$data,nrow)))
    loglik <- logLik(object,...)

    df <- attributes(loglik)$df
    nobs <- attributes(loglik)$nall*length(endogenous(object))
    myAIC <- -2*(loglik - df); attributes(myAIC) <- NULL
    myBIC <- -2*loglik + df*log(nobs); attributes(myBIC) <- NULL

    xconstrain <- intersect(unlist(lapply(constrain(object),function(z) attributes(z)$args)),manifest(object))

    l2D <- mean(object$opt$grad^2)
    S <- vcov(object)
    rnkV <- tryCatch(qr(S)$rank,error=function(...) 0)
    minSV <- attr(S,"minSV")
    condnum <- tryCatch(condition(vcov(object)),error=function(...) NULL)

    if (((object$estimator=="gaussian" & class(object)[1]!="lvm.missing") | chisq) & length(xconstrain)==0 ) {
        res <- list(fit=compare(object), n=n, logLik=loglik, BIC=myBIC, AIC=myAIC)
        q <- res$fit$statistic
        qdf <- res$fit$parameter
        if (all) {
            m0 <- lvm(manifest(object)); exogenous(m0) <- NULL
            e0 <- estimate(m0,model.frame(object))
            g0 <- gof(e0)
            qbaseline <- g0$fit$statistic
            qdfbaseline <- g0$fit$parameter
            CFI <- ((qbaseline-qdfbaseline) - (q-qdf))/(qbaseline-qdfbaseline)
            NFI <- (qbaseline-q)/qbaseline
            TLI <- (qbaseline/qdfbaseline-q/qdf)/(qbaseline/qdfbaseline-1)
            S <- object$data$S
            mu <- object$data$mu
            C <- modelVar(object)$C
            xi <- as.vector(modelVar(object)$xi)
            if (is.null(S)) S <- cov(model.frame(object))
            if (is.null(mu)) mu <- colMeans(model.frame(object))
            L <- diag(S)^0.5
            R <- (diag(1/L))%*%(S-C)%*%(diag(1/L))
            R2 <- (mu-xi)/L
            SRMR <- mean(c(R[upper.tri(R,diag=TRUE)],R2)^2)^0.5
            res <- c(res,list(CFI=CFI,NFI=NFI,TLI=TLI,C=C,S=S,SRMR=SRMR))
            ## idx <- index(object)$endo.idx
            ## if (length(latent(object))>0) {
            ##   SRMR.endo <- mean(c(R[idx,idx][upper.tri(R[idx,idx],diag=TRUE)],R2[idx])^2)^0.5
            ##   res <- c(res,list("SRMR(endogenous)"=SRMR.endo))
            ## }
        }
        if (rnkV==ncol(vcov(object)) && (!is.null(minSV) && minSV>1e-12)) {

            rmseafun <- function(...) {
                epsilon <- function(lambda) sapply(lambda,function(x)
                    ifelse(x>0 & qdf>0,sqrt(x/(qdf*(n))),0)) ## n-1,n vs. n-df
                opf <- function(l,p) suppressWarnings(p-pchisq(q,df=qdf,ncp=l))
                alpha <- (1-level)/2
                RMSEA <- epsilon(q-qdf)
                B <- max(q-qdf,0)
                lo <- hi <- list(root=0)
                if (RMSEA>0 && opf(0,p=1-alpha)<0) {
                    hi <- uniroot(function(x) opf(x,p=1-alpha),c(0,B))
                }
                if (opf(B,p=alpha)<0) {
                    lo <- uniroot(function(x) opf(x,p=alpha),c(B,n))
                }
                ci <- c(epsilon(c(hi$root,lo$root)))
                RMSEA <- c(RMSEA=RMSEA,ci);
                names(RMSEA) <- c("RMSEA",paste0(100*c(alpha,(1-alpha)),"%"))
                pval <- pchisq(q,qdf,(n*qdf*rmsea.threshold^2),lower.tail=FALSE)
                res <- list(aa=((q-qdf)/(2*qdf)^0.5),RMSEA=RMSEA, level=level, rmsea.threshold=rmsea.threshold, pval.rmsea=pval)
                return(res)
            }
            rmseaval <- tryCatch(rmseafun(),error=function(e) NULL)
            res <- c(res,rmseaval)
        }
    } else {
        res <- list(n=n, logLik=loglik, BIC=myBIC, AIC=myAIC)
    }

    res <- c(res, L2score=l2D, rankV=rnkV, cond=condnum, k=nrow(vcov(object)))
    class(res) <- "gof.lvmfit"
    return(res)
}

##' @export
print.gof.lvmfit <- function(x,optim=TRUE,...) {
    if (!is.null(x$n)) {
        with(x,
             cat("\n Number of observations =", n, "\n"))
    }
    if (is.null(x$fit)) {
        with(x,
             cat(" Log-Likelihood =", logLik, "\n"))
    }
    with(x,  cat(" BIC =", BIC, "\n",
                 "AIC =", AIC, "\n"))
    if (!is.null(x$fit))
        with(x,
             cat(" log-Likelihood of model =", fit$estimate[1], "\n\n",
                 "log-Likelihood of saturated model =", fit$estimate[2], "\n",
                 "Chi-squared statistic: q =", fit$statistic,
                 ", df =", fit$parameter,
                 "\n  P(Q>q) =", fit$p.value, "\n"))
    if (!is.null(x$RMSEA)) {
        rr <- round(x$RMSEA*10000)/10000
        rmsea <- paste0(rr[1]," (",rr[2],";",rr[3],")")
        cat("\n RMSEA (",x$level*100,"% CI): ", rmsea,"\n",sep="")
        cat("  P(RMSEA<",x$rmsea.threshold,")=",  x$pval.rmsea,"\n",sep="")
    }
    for (i in c("TLI","CFI","NFI","SRMR","SRMR(endogenous)"))
        if (!is.null(x[[i]])) cat("", i,"=",x[[i]],"\n")

    if (optim) {
        cat("\nrank(Information) = ",x$rankV," (p=", x$k,")\n",sep="")
        cat("condition(Information) = ",x$cond,"\n",sep="")
        cat("mean(score^2) =",x$L2score,"\n")
    }

    invisible(x)
}



## gof.multigroupfit <- function(object,...) {
##   L0 <- logLik(object); df0 <- attributes(L0)$df
##   L1 <- logLik(object,type="sat"); df1 <- attributes(L1)$df

##   df <- df1-df0; names(df) <- "df"
##   Q <- -2*(L0-L1); attributes(Q) <- NULL; names(Q) <- "chisq";
##   pQ <- pchisq(Q,df,lower.tail=FALSE)
##   values <- c(L0,L1); names(values) <- c("log likelihood (model)", "log likelihood (saturated model)")
##   res <- list(statistic = Q, parameter = df,
##               p.value=pQ, method = "Likelihood ratio test",
##               estimate = values)
##   class(res) <- "htest"
##   return(res)
## }