.packageName <- "sem"
## last modified 9 November 02 by J. Fox

mod.indices<-function(model, ...){
    UseMethod("mod.indices")
    }

mod.indices.sem <- function(model, ...){        
    accumulate <- function(A, B, C, D, d) {
        res <- matrix(0, d^2, d^2)    
        for (ii in 1:(d^2)){
            for (jj in ii:(d^2)){
                g <- 1 + (ii - 1) %% d
                h <- 1 + (ii - 1) %/% d
                i <- 1 + (jj - 1) %% d
                j <- 1 + (jj - 1) %/% d
                res[ii, jj] <- res[jj, ii] <- A[g, i] * B[h, j] + C[g, j] * D[h, i]
                }
            }
        res
        }
    accumulate.asym <- function(A, B, C, D, d) {
        res <- matrix(0, d^2, d^2)    
        for (ii in 1:(d^2)){
            for (jj in 1:(d^2)){
                g <- 1 + (ii - 1) %% d
                h <- 1 + (ii - 1) %/% d
                i <- 1 + (jj - 1) %% d
                j <- 1 + (jj - 1) %/% d
                res[ii, jj] <- A[g, i] * B[h, j] + C[g, j] * D[h, i]
                }
            }
        res
        }        
    A <- model$A
    P <- model$P
    S <- model$S
    C <- model$C
    J <- model$J
    m <- model$m
    NM <- model$N - 1
    vars <- model$var.names    
    I.Ainv <- solve(diag(m) - A) 
    Cinv <- solve(C)    
    AA <- t(I.Ainv) %*% t(J)
    BB <- J %*% I.Ainv %*% P %*% t(I.Ainv)
    grad.P <- AA %*% Cinv %*% (C - S) %*% Cinv %*% t(AA)
    grad.A <- AA %*% Cinv %*% (C - S) %*% Cinv %*% BB 
    grad <- c(grad.A, grad.P)        
    dF.dBdB <- accumulate(AA %*% Cinv %*% t(AA), t(BB) %*% Cinv %*% BB,
                    AA %*% Cinv %*% BB, t(BB) %*% Cinv %*% t(AA), m)                
    dF.dPdP <- accumulate(AA %*% Cinv %*% t(AA), AA %*% Cinv %*% t(AA),
                    AA %*% Cinv %*% t(AA), AA %*% Cinv %*% t(AA), m)                
    dF.dBdP <- accumulate.asym(AA %*% Cinv %*% t(AA), t(BB) %*% Cinv %*% t(AA),
                    AA %*% Cinv %*% t(AA), t(BB) %*% Cinv %*% t(AA), m)    
    ram <- model$ram   
    fixed <- ram[,4] == 0
    one.head <- ram[,1] == 1
    one.free <- which( (!fixed) & one.head )
    two.free <- which( (!fixed) & (!one.head) )
    arrows.1.free <- ram[one.free,c(2,3)]
    arrows.2.free <- ram[two.free,c(2,3)]
    posn.matrix <- matrix(1:(m^2), m, m)
    posn.free <- c(posn.matrix[arrows.1.free], 
                (m^2) + posn.matrix[arrows.2.free])                    
    Hessian <- rbind( cbind(dF.dBdB,    dF.dBdP),
                      cbind(t(dF.dBdP), dF.dPdP))
    hessian <- Hessian[posn.free, posn.free]
    E.inv <- (2/NM) * solve(hessian)                      
    par.change <- mod.indices <- rep(0, 2*(m^2))                
    for (i in 1:(2*(m^2))) {
        k <- Hessian[i, i]
        d <- Hessian[i, posn.free]
     #   mod.indices[i] <- NM * grad[i]^2 / (k - d %*% E.inv %*% d)
        par.change[i] <- -grad[i]/ (k - d %*% E.inv %*% d)
        mod.indices[i] <- -NM * grad[i] * par.change[i]
        }
    mod.A <- matrix(mod.indices[1:(m^2)], m, m)
    mod.P <- matrix(mod.indices[-(1:(m^2))], m, m)
    par.A <- matrix(par.change[1:(m^2)], m, m)
    par.P <- matrix(par.change[-(1:(m^2))], m, m)
    mod.A[arrows.1.free] <- NA
    par.A[arrows.1.free] <- NA
    mod.P[arrows.2.free] <- NA
    par.P[arrows.2.free] <- NA
    rownames(mod.A) <- colnames(mod.A) <- vars
    rownames(mod.P) <- colnames(mod.P) <- vars
    rownames(par.A) <- colnames(par.A) <- vars
    rownames(par.P) <- colnames(par.P) <- vars
    result <- list(mod.A=mod.A, mod.P=mod.P, par.A=par.A, par.P=par.P)
    class(result) <- "sem.modind"
    result
    }

summary.sem.modind <- function(object, round=2, 
    print.matrices=c('both', 'par.change', 'mod.indices'), ...) {
    print.matrices <- match.arg(print.matrices)
    if (print.matrices != "mod.indices"){ 
        cat("\n Parameter change: A matrix\n")
        print(object$par.A)
        }
    if (print.matrices != "par.change"){ 
        cat("\n Modification indices: A matrix\n")
        print(round(object$mod.A, round))
        }
    if (print.matrices != "mod.indices"){ 
        cat("\n Parameter change: P matrix\n")
        print(object$par.P)
        }
    if (print.matrices != "par.change"){ 
        cat("\n Modification indices: P matrix\n")
        print(round(object$mod.P, round))
        }
    invisible(NULL)
    }

print.sem.modind <- function(x, n.largest=5, ...){
    cat("\n", n.largest, "largest modification indices, A matrix:\n")
    mod.A <- as.vector(x$mod.A)
    names <- rownames(x$mod.A)
    names(mod.A) <- outer(names, names, paste, sep=":")
    print(rev(sort(mod.A))[1:n.largest])
    cat("\n ", n.largest, "largest modification indices, P matrix:\n")
    mod.P <- as.vector(x$mod.P)
    names(mod.P) <- outer(names, names, paste, sep=":")
    print(rev(sort(mod.P[lower.tri(x$mod.P, diag=TRUE)]))[1:n.largest])
    }
# last modified 9 November 02 by J. Fox

path.diagram <- function(model, ...){
    UseMethod("path.diagram")
    }

path.diagram.sem <- function(model, out.file, min.rank=NULL, max.rank=NULL, same.rank=NULL,
    variables=model$var.names, parameters=rownames(object$ram), ignore.double=TRUE,
    edge.labels=c("names", "values"), size=c(8,8), node.font=c("Helvetica", 14),
    edge.font=c("Helvetica", 10), rank.direction=c("LR", "TB"), ...){
    if(!missing(out.file)){
        handle <- file(out.file, "w")
        on.exit(close(handle))
        }
        else handle <- stdout()
    edge.labels <- match.arg(edge.labels)
    rank.direction <- match.arg(rank.direction)
    cat(file=handle, paste('digraph "', deparse(substitute(model)),
        '" {\n', sep=""))
    cat(file=handle, paste('  rankdir=', rank.direction, ';\n', sep=""))
    cat(file=handle, paste('  size="',size[1],',',size[2],'";\n', sep=""))
    cat(file=handle, paste('  node [fontname="', node.font[1], 
        '" fontsize=', node.font[2], ' shape=box];\n', sep=""))
    cat(file=handle, paste('  edge [fontname="', edge.font[1],
        '" fontsize=', edge.font[2], '];\n', sep=""))
    cat(file=handle, '  center=1;\n')
    if (!is.null(min.rank)){
        min.rank <- paste('"', min.rank, '"', sep="")
        min.rank <- gsub(',', '" "', gsub(' ', '', min.rank))
        cat(file=handle, paste('  {rank=min ', min.rank, '}\n', sep=""))
        }
    if (!is.null(max.rank)){
        max.rank <- paste('"', max.rank, '"', sep="")
        max.rank <- gsub(',', '" "', gsub(' ', '', max.rank))
        cat(file=handle, paste('  {rank=max ', max.rank,'}\n', sep=""))
        }
    if (!is.null(same.rank)){
        for (s in 1:length(same.rank)){
            same <- paste('"', same.rank[s], '"', sep="")
            same <- gsub(',', '" "', gsub(' ', '', same))
            cat(file=handle, paste('  {rank=same ', same,'}\n', sep=""))
            }
        }
    latent <- variables[-(1:model$n)]
    for (lat in latent){
        cat(file=handle, paste('  "', lat, '" [shape=ellipse]\n', sep=""))
        }
    ram <- model$ram
    ram[names(model$coeff), 5] <- model$coeff
    rownames(ram)[model$fixed] <- ram[model$fixed,5]
    names <- rownames(ram)
    values <- ram[,5]
    heads <- ram[,1]
    to <- ram[,2]
    from <- ram[,3]
    labels <- if (edge.labels == "names") names else values
    direction <- ifelse((heads == 2), ' dir=both', "")
    for (par in 1:nrow(ram)){
        if((!ignore.double) || (heads[par] == 1))
        cat(file=handle, paste('  "', 
            variables[from[par]], '" -> "', variables[to[par]], 
            '" [label="', labels[par], '"', direction[par], '];\n', sep=""))
        }
    cat(file=handle, '}\n')
    }
# last modified 27 July 2002 by J. Fox

print.sem <- function(x, ...) {
    n <- x$n
    t <- x$t
    n.fix <- x$n.fix
    df <- n*(n + 1)/2 - t - n.fix*(n.fix + 1)/2
    cat("\n Model Chisquare = ", x$criterion * (x$N - 1), 
        "  Df = ", df, "\n\n")
    print(x$coeff)
    cat("\n Iterations = ", x$iterations, "\n")
    if (!is.null(x$aliased)) cat("\n Aliased parameters:", x$aliased, "\n")
    invisible(x)
    }
# last modified 27 March 02 by J. Fox

ram<-function(object, digits=5, startvalues=FALSE){
    old.digits <- options(digits=digits)
    on.exit(options(old.digits))
    var.names <- rownames(object$A)
    ram <- object$ram
    if (!startvalues) colnames(ram) <- c(colnames(ram)[1:4], "estimate")
    par <- object$coeff
    par.names <- rep(" ", nrow(ram))
    t <- object$t
    for (i in 1:t) {
        which.par <- ram[,4] == i
        if (!startvalues)  ram[which.par, 5] <- par[i]
        par.names[which.par] <- names(par)[i]
        }
    par.code <- paste(var.names[ram[,2]], c('<---', '<-->')[ram[,1]],
                    var.names[ram[,3]])
    ram <- data.frame(ram, arrow = par.code)
    print(ram, rowlab=par.names)
    }
# last modified 14 April 2001 by J. Fox

residuals.sem <- function(object, ...){
    object$S - object$C
    }    

standardized.residuals <- function(object, ...){
    UseMethod("standardized.residuals")
    }

standardized.residuals.sem <- function(object, ...){
    res <- residuals(object)
    s <- diag(object$S)
    res/sqrt(outer(s, s))
    }

normalized.residuals <- function(object, ...){
    UseMethod("normalized.residuals")
    }
    
normalized.residuals.sem <- function(object, ...){
    res <- residuals(object)
    N <- object$N
    C <- object$C
    c <- diag(C)
    res/sqrt((outer(c,c) + C^2)/N)
    }
# last modified 7 Mar 04 by J. Fox

sem <- function(ram, ...){
    if (is.character(ram)) class(ram) <- 'mod'
    UseMethod('sem', ram)
    }

sem.mod <- function (ram, S, N, obs.variables=rownames(S), fixed.x=NULL, debug=FALSE, ...){
    parse.path <- function(path) {
        path.1 <- gsub('-', '', gsub(' ','', path))
        direction <- if (regexpr('<>', path.1) > 0) 2 
            else if (regexpr('<', path.1) > 0) -1
            else if (regexpr('>', path.1) > 0) 1
            else stop(paste('ill-formed path:', path))
        path.1 <- strsplit(path.1, '[<>]')[[1]]
        list(first=path.1[1], second=path.1[length(path.1)], direction=direction)
        }
    if ((!is.matrix(ram)) | ncol(ram) != 3) stop ('ram argument must be a 3-column matrix')
    startvalues <- as.numeric(ram[,3])
    par.names <- ram[,2]
    n.paths <- length(par.names)
    heads <- from <- to <- rep(0, n.paths)
    for (p in 1:n.paths){
        path <- parse.path(ram[p,1])
        heads[p] <- abs(path$direction)
        to[p] <- path$second
        from[p] <- path$first
        if (path$direction == -1) {
            to[p] <- path$first
            from[p] <- path$second
            }
        }
    ram <- matrix(0, p, 5)
    all.vars <- unique(c(to, from))
    latent.vars <- setdiff(all.vars, obs.variables)
    vars <- c(obs.variables, latent.vars)
    pars <- na.omit(unique(par.names))
    ram[,1] <- heads
    ram[,2] <- apply(outer(vars, to, '=='), 2, which)
    ram[,3] <- apply(outer(vars, from, '=='), 2, which)   
    par.nos <- apply(outer(pars, par.names, '=='), 2, which)
    ram[,4] <- unlist(lapply(par.nos, function(x) if (length(x) == 0) 0 else x))
    ram[,5]<- startvalues
    colnames(ram) <- c('heads', 'to', 'from', 'parameter', 'start')
    if (!is.null(fixed.x)) fixed.x <- apply(outer(vars, fixed.x, '=='), 2, which)
    n <- length(obs.variables)
    m <- length(all.vars)
    t <- length(pars)
    if (debug) {
        cat('\n observed variables:\n') 
        print(paste(paste(1:n,':', sep=''), obs.variables, sep=''))
        cat('\n')
        if (m > n){ 
            cat('\n latent variables:\n')
            print(paste(paste((n+1):m,':', sep=''), latent.vars, sep=''))
            cat('\n')
            }
        cat('\n parameters:\n') 
        print(paste(paste(1:t,':', sep=''), pars, sep=''))
        cat('\n\n RAM:\n')
        print(ram)
        }
    sem(ram=ram, S=S, N=N, param.names=pars, var.names=vars, fixed.x=fixed.x,
        debug=debug, ...)
    }
     

sem.default <- function(ram, S, N, param.names=paste('Param', 1:t, sep=''), 
    var.names=paste('V', 1:m, sep=''), fixed.x=NULL, debug=FALSE,
    analytic.gradient=TRUE, warn=FALSE, maxiter=500, par.size=c('ones', 'startvalues'), 
    refit=TRUE, start.tol=1E-6, ...){
    ord <- function(x) 1 + apply(outer(unique(x), x, "<"), 2, sum)
    is.triangular <- function(X) {
        is.matrix(X) && (nrow(X) == ncol(X)) && 
            (all(0 == X[upper.tri(X)])) || (all(0 == X[lower.tri(X)]))
        }    
    is.symmetric <- function(X) {
        is.matrix(X) && (nrow(X) == ncol(X)) && all(X == t(X))
        }
    if (is.triangular(S)) S <- S + t(S) - diag(diag(S))
    if (!is.symmetric(S)) stop('S must be a square triangular or symmetric matrix')
    if ((!is.matrix(ram)) | ncol(ram) != 5 | (!is.numeric(ram)))
        stop ('ram argument must be a 5-column numeric matrix')
    par.size <- if (missing(par.size)) {
        range <- range(diag(S))
        if (range[2]/range[1] > 100) 'startvalues' else 'ones'
        }
        else match.arg(par.size)
    n <- nrow(S)
    observed <- 1:n
    n.fix <- length(fixed.x)
    if (!is.null(fixed.x)){
        for (i in 1:n.fix){
            for (j in 1:i){
                ram <- rbind(ram, c(2, fixed.x[i], fixed.x[j], 
                    0, S[fixed.x[i], fixed.x[j]]))
                }
            }
        }
    m <- max(ram[,2])
    missing.variances <- setdiff(1:m, ram[,2][ram[,2] == ram[,3]])
    if (length(missing.variances) > 0) warning(paste("The following variables have no variance or error-variance parameter (double-headed arrow):\n",
        paste(var.names[missing.variances], collapse=", "), 
        "\nThe model is almost surely misspecified; check also for missing covariances.\n"))
    t <- max(ram[,4])
    df <- n*(n + 1)/2 - t - n.fix*(n.fix + 1)/2
    if (df < 0) stop(paste("The model has negative degrees of freedom =", df))
    J <- matrix(0, n, m)
    correct <- matrix(2, m, m)
    diag(correct) <- 1
    J[cbind(1:n, observed)]<-1
    par.posn <-  sapply(1:t, function(i) which(ram[,4] == i)[1])
    colnames(ram)<-c("heads", "to", "from", "parameter", "start value")
    rownames(ram)<-rep("",nrow(ram))
    rownames(ram)[par.posn]<-param.names
    fixed <- ram[,4] == 0
    sel.free <- ram[,4]
    sel.free[fixed] <- 1
    one.head <- ram[,1] == 1
    one.free <- which( (!fixed) & one.head )
    two.free <- which( (!fixed) & (!one.head) )
    arrows.1 <- ram[one.head, c(2,3), drop=FALSE]
    arrows.2 <- ram[!one.head, c(2,3), drop=FALSE]
    arrows.2t <- ram[!one.head, c(3,2), drop=FALSE]
    arrows.1.free <- ram[one.free,c(2,3), drop=FALSE]
    arrows.2.free <- ram[two.free,c(2,3), drop=FALSE]
    sel.free.1 <- sel.free[one.free]
    sel.free.2 <- sel.free[two.free]
    unique.free.1 <- unique(sel.free.1)
    unique.free.2 <- unique(sel.free.2)
    start <- if (any(is.na(ram[,5][par.posn]))) startvalues(S, ram, debug=debug, tol=start.tol)
        else ram[,5][par.posn]
    objective.1 <- function(par){
        A <- P <- matrix(0, m, m)
        val <- ifelse (fixed, ram[,5], par[sel.free])
        A[arrows.1] <- val[one.head]
        P[arrows.2t] <- P[arrows.2] <- val[!one.head]
        I.Ainv <- solve(diag(m) - A)
        C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)
        Cinv <- solve(C)
        f <- sum(diag(S %*% Cinv)) + log(det(C))
        f
        }
    objective.2 <- function(par){
        A <- P <- matrix(0, m, m)
        val <- ifelse (fixed, ram[,5], par[sel.free])
        A[arrows.1] <- val[one.head]
        P[arrows.2t] <- P[arrows.2] <- val[!one.head]
        I.Ainv <- solve(diag(m) - A)
        C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)
        Cinv <- solve(C)
        f <- sum(diag(S %*% Cinv)) + log(det(C))
        grad.P <- correct * t(I.Ainv) %*% t(J) %*% Cinv %*% (C - S) %*% Cinv %*% J %*% I.Ainv
        grad.A <- grad.P %*% P %*% t(I.Ainv)        
        gradient <- rep(0, t)
        gradient[unique.free.1] <- tapply(grad.A[arrows.1.free],ram[ram[,1]==1 & ram[,4]!=0, 4], sum)
        gradient[unique.free.2] <- tapply(grad.P[arrows.2.free],ram[ram[,1]==2 & ram[,4]!=0, 4], sum)
        attributes(f) <- list(C=C, A=A, P=P, gradient=gradient)
        f
        }
    if (!warn){
        save.warn <- options(warn=-1)
        on.exit(options(save.warn))
        }
    typsize <- if (par.size == 'startvalues') abs(start) else rep(1,t)
    res <- nlm(if (analytic.gradient) objective.2 else objective.1, 
        start, hessian=TRUE, iterlim=maxiter, print.level=if(debug) 2 else 0,
        typsize=typsize, ...)
    convergence <- res$code
    if (!warn) options(save.warn)
    par <- res$estimate
    names(par) <- param.names
    result <- list()
    result$var.names <- var.names
    obj <- objective.2(par)
    ram[par.posn, 5] <- start
    par.code <- paste(var.names[ram[,2]], c('<---', '<-->')[ram[,1]],
    var.names[ram[,3]])
    result$ram <- ram
    result$coeff <- par
    result$criterion <-  c(obj) - n - log(det(S))
    rownames(S) <- colnames(S) <- var.names[observed]
    result$S <- S
    result$J <- J
    C <- attr(obj, "C")
    rownames(C) <- colnames(C) <- var.names[observed]
    result$C <- C
    A <- attr(obj, "A")
    rownames(A) <- colnames(A) <- var.names
    result$A <- A
    P <- attr(obj, "P")
    rownames(P) <- colnames(P) <- var.names
    result$P <- P
    result$n.fix <- n.fix
    result$n <- n
    result$N <- N
    result$m <- m
    result$t <- t
    result$par.posn <- par.posn
    result$convergence <- convergence
    result$iterations <- res$iterations
    if (convergence > 2) 
        warning(paste('Optimization may not have converged; nlm return code = ',
            res$code, '. Consult ?nlm.\n', sep=""))
    qr.hess <- try(qr(res$hessian), silent=TRUE)
    if (class(qr.hess) == "try-error")
        warning("Could not compute QR decomposition of Hessian.\nOptimization probably did not converge.\n")
    else if (qr.hess$rank < t){ 
        warning(' singular Hessian: model is probably underidentified.\n')
        cov <- matrix(NA, t, t)
        colnames(cov) <- rownames(cov) <- param.names
        result$cov <- cov
        which.aliased <- qr.hess$pivot[-(1:qr.hess$rank)]
        aliased <- param.names[which.aliased]
        position.aliased <- is.element(ram[,4], which.aliased)
        if (refit){
            warning(' refitting without aliased parameters.\n')
            ram.refit <- ram[!position.aliased,]
            par.refit <- ram.refit[,4]
            par.refit[par.refit != 0] <- ord(par.refit[par.refit != 0])
            ram.refit[,4] <- par.refit
            iterations <- result$iterations
            result <- Recall(ram.refit, S, N, 
                param.names=param.names[-which.aliased], var.names=var.names,
                debug=debug, analytic.gradient=analytic.gradient, 
                warn=warn, maxiter=maxiter, par.size=par.size, refit=TRUE, ...)
            result$iterations <- iterations + result$iterations
            aliased <- c(aliased, result$aliased)
            }
        result$aliased <- aliased
        }
    else {
        cov <- (2/(N - 1)) * solve(res$hessian)
        colnames(cov) <- rownames(cov) <- param.names
        result$cov <- cov
        if (any(diag(cov) < 0)) {
            result$aliased <- c(param.names[diag(cov) < 0], result$aliased)
            warning("Negative parameter variances.\nModel is probably underidentified.\n")
            }
        }
    class(result) <- "sem"
    result
    }
# last modified 25 July 2001 by J. Fox

standardized.coefficients <- function(object, digits=5){
    old.digits <- options(digits=digits)
    on.exit(options(old.digits))
    P <- object$P
    A <- object$A
    t <- object$t
    par <- object$coeff
    par.posn <- object$par.posn
    IAinv <- solve(diag(nrow(A)) - A)
    C <- IAinv %*% P %*% t(IAinv)
    ram <- object$ram
    par.names <- rep(' ', nrow(ram))
    for (i in 1:t) {
        which.par <- ram[,4] == i
        ram[which.par, 5] <- par[i]
        par.names[which.par] <- names(par)[i]
        }
    one.head <- ram[,1] == 1
    coeff <- ram[one.head, 5]
    coeff <- coeff *
        sqrt(diag(C[ram[one.head, 3], ram[one.head, 3]])/diag(C[ram[one.head, 2], ram[one.head, 2]]))
    var.names <- rownames(A)   
    par.code <- paste(var.names[ram[one.head,2]], c('<---', '<-->')[ram[one.head,1]],
                    var.names[ram[one.head,3]])
    coeff <- data.frame(par.names[one.head], coeff, par.code)
    names(coeff) <- c(" ", "Std. Estimate", " ")
    print(coeff, rowlab=rep(" ", nrow(coeff)))
    }

std.coef <- function (...){
    standardized.coefficients(...)
    }
# last modified 28 Feb 2004 by J. Fox

startvalues <- function(S, ram, debug=FALSE, tol=1E-6){
    n <- nrow(S) 
    observed <- 1:n       
    m <- max(ram[,2])            
    t <- max(ram[,4])   
    s <- sqrt(diag(S))
    R <- S/outer(s,s)
    latent<-(1:m)[-observed]
    par.posn <-  sapply(1:t, function(i) which(ram[,4] == i)[1])    
    one.head <- ram[,1] == 1
    start <- (ram[,5])[par.posn]
    A.pat <-matrix(FALSE, m, m)
    A.pat[ram[one.head, c(2,3), drop=FALSE]] <- TRUE
    P.pat <- C <- matrix(0, m, m)
    P.pat[ram[!one.head, c(2,3), drop=FALSE]] <- P.pat[ram[!one.head, c(3,2), drop=FALSE]] <- 1
    C[observed, observed] <- R
    for (l in latent) {
        indicators <- A.pat[observed, l]
        for (j in observed){
            C[j, l] <- C[l, j] <- if (!any(indicators)) runif(1, .3, .5)
                else {   
                        numerator <- sum(R[j, observed[indicators]])
                        denominator <- sqrt(sum(R[observed[indicators], observed[indicators]]))
                        numerator/denominator
                    }
            }
        }
    for (l in latent){
        for (k in latent){
            C[l, k] <- C[k,l] <- if (l==k) 1 else {
                                indicators.l <- A.pat[observed, l]
                                indicators.k <- A.pat[observed, k]
                                if ((!any(indicators.l)) | (!any(indicators.k))) runif(1, .3, .5) else {
                                    numerator <- sum(R[observed[indicators.l], observed[indicators.k]])
                                    denominator <- sqrt( sum(R[observed[indicators.l], observed[indicators.l]])
                                        * sum(R[observed[indicators.k], observed[indicators.k]]))
                                    numerator/denominator}
                                    }
            }
        }
    A <- matrix(0, m, m)
    for (j in 1:m){
        ind <- A.pat[j,]
        if (!any(ind)) next
        A[j, ind] <- solve(C[ind, ind]) %*% C[ind, j]
        }
    A[observed,] <- A[observed,]*matrix(s, n, m)
    A[,observed] <- A[,observed]/matrix(s, m, n, byrow=TRUE)
    C[observed,] <- C[observed,]*matrix(s, n, m)
    C[,observed] <- C[,observed]*matrix(s, m, n, byrow=TRUE)
    P <- (diag(m) - A) %*% C %*% t(diag(m) - A)
    P <- P.pat * P
    for (par in 1:t){
        if (!is.na(start[par])) next
        posn <- par.posn[par]
        if (ram[posn, 1] == 1) start[par] <- A[ram[posn, 2], ram[posn, 3]]
            else start[par] <- P[ram[posn, 2], ram[posn, 3]]
        if (abs(start[par]) < tol) start[par] <- tol
        }
    if (debug){
        cat('\nStart values:\n')
        print(start)
        cat('\n')
        }
    start
    }
# last modified 28 Feb 2004 by J. Fox

summary.sem <- function(object, digits=5, conf.level=.90, ...) {
    norm.res <- normalized.residuals(object)
    se <- sqrt(diag(object$cov))
    z <- object$coeff/se
    n.fix <- object$n.fix
    n <- object$n
    t <- object$t
    S <- object$S
    C <- object$C
    N <- object$N
    df <- n*(n + 1)/2 - t - n.fix*(n.fix + 1)/2
    invC <- solve(C)
    CSC <- invC %*% (S - C)
    CSC <- CSC %*% CSC
    CS <- invC %*% S
    CS <- CS %*% CS
    GFI <- 1 - sum(diag(CSC))/sum(diag(CS))
    AGFI <- if (df > 0) 1 - (n*(n + 1)/(2*df))*(1 - GFI)
     else NA
    chisq <- object$criterion * (N - 1)
    RMSEA <- sqrt(max(object$criterion/df - 1/(N - 1), 0))
    tail <- (1 - conf.level)/2 
    max <- 1000
    while (max > 1){
        res <- optimize(function(lam) (tail - pchisq(chisq, df, ncp=lam))^2, interval=c(0, max))
        if (sqrt(res$objective) < tail/100) break
        max <- max/2
        }
    lam.U <- if (max <= 1) NA else res$minimum
    while (max > 1){
        res <- optimize(function(lam) (1 - tail - pchisq(chisq, df, ncp=lam))^2, interval=c(0, max))
        if (sqrt(res$objective) < tail/100) break
        max <- max/2
        }
    lam.L <- if (max <= 1) NA else res$minimum
    RMSEA.U <- sqrt(lam.U/((N - 1)*df))
    RMSEA.L <- sqrt(lam.L/((N - 1)*df))
    RMSEA <- c(RMSEA, RMSEA.L, RMSEA.U, conf.level)
    var.names <- rownames(object$A)
    ram <- object$ram[object$par.posn,]
    par.code <- paste(var.names[ram[,2]], c('<---', '<-->')[ram[,1]],
                    var.names[ram[,3]])
    coeff <- data.frame(object$coeff, se, z, 2*(1 - pnorm(abs(z))), par.code)
    names(coeff) <- c("Estimate", "Std Error", "z value", "Pr(>|z|)", " ")
    row.names(coeff) <- names(object$coeff)
    BIC <- if (df > 0) chisq - df * log(N*n) else NA
    ans <- list(chisq=chisq, df=df, GFI=GFI, AGFI=AGFI, RMSEA=RMSEA, BIC=BIC, 
        norm.res=norm.res, coeff=coeff, digits=digits, 
        iterations=object$iterations, aliased=object$aliased)
    class(ans) <- "summary.sem"
    ans
    }
    
print.summary.sem <- function(x, ...){
    old.digits <- options(digits=x$digits)
    on.exit(options(old.digits))
    cat("\n Model Chisquare = ", x$chisq, "  Df = ", x$df, 
        "Pr(>Chisq) =", if (x$df > 0) 1 - pchisq(x$chisq, x$df)
            else NA)
    cat("\n Goodness-of-fit index = ", x$GFI)
    cat("\n Adjusted goodness-of-fit index = ", x$AGFI)
    cat("\n RMSEA index =  ", x$RMSEA[1],
        "   ", 100*x$RMSEA[4], " \% CI: (", x$RMSEA[2], ", ", x$RMSEA[3],")", sep="")
    cat("\n BIC = ", x$BIC, "\n")
    cat("\n Normalized Residuals\n")
    print(summary(as.vector(x$norm.res)))
    cat("\n Parameter Estimates\n")
    print(x$coeff)
    cat("\n Iterations = ", x$iterations, "\n")
    if (!is.null(x$aliased)) cat("\n Aliased parameters:", x$aliased, "\n")
    invisible(x)
    }
# Two-Stage Least Squares
#   John Fox

# last modified 17 Nov 03 by J. Fox

tsls <- function(y, ...){
    UseMethod("tsls")
    }

tsls.default <- function (y, X, Z, names=NULL, ...) {
    n <- length(y)
    p <- ncol(X)
    invZtZ <- solve(crossprod(Z))
    XtZ <- crossprod(X, Z)
    V <- solve(XtZ %*% invZtZ %*% t(XtZ))
    b <- V %*% XtZ %*% invZtZ %*% crossprod(Z, y)
    residuals <- y - X %*% b
    s2 <- sum(residuals^2)/(n - p)
    V <- s2*V
    result<-list()
    result$n <- n
    result$p <- p
    b <- as.vector(b)
    names(b) <- names
    result$coefficients <- b
    rownames(V) <- colnames(V) <- names
    result$V <- V
    result$s <- sqrt(s2)
    result$residuals <- as.vector(residuals)
    result$response <- y
    result$model.matrix <- X
    result$instruments <- Z
    result
    }

    
tsls.formula <- function (formula, instruments, data, subset, na.action, contrasts = NULL, ...) {
    if (missing(na.action)) 
        na.action <- options()$na.action
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval(m$data, sys.frame(sys.parent())))) 
        m$data <- as.data.frame(data)
    response.name <- deparse(formula[[2]])
    form <- as.formula(paste(response.name, "~", deparse(formula[[3]]), 
        "+", deparse(instruments[[2]])))
    m$formula <- form
    m$instruments <- m$contrasts <- NULL
    m[[1]] <- as.name("model.frame")
    mf <- eval(m, sys.frame(sys.parent()))
    na.act <- attr(mf, "na.action")
    Z <- model.matrix(instruments, data = mf, contrasts)
    y <- mf[, response.name]
    X <- model.matrix(formula, data = mf, contrasts)
    result <- tsls(y, X, Z, colnames(X))
    result$response.name <- response.name
    result$formula <- formula
    result$instruments <- instruments
    if (!is.null(na.act)) 
        result$na.action <- na.act
    class(result) <- "tsls"
    result
    }


print.tsls <- function(x, ...){
    cat("\nModel Formula: ")
    print(x$formula)
    cat("\nInstruments: ")
    print(x$instruments)
    cat("\nCoefficients:\n")
    print(x$coefficients)
    cat("\n")
    invisible(x)
    }
    
    
summary.tsls <- function(object, digits=4, ...){
    save.digits <- unlist(options(digits=digits))
    on.exit(options(digits=save.digits))
    cat("\n 2SLS Estimates\n")
    cat("\nModel Formula: ")
    print(object$formula)
    cat("\nInstruments: ")
    print(object$instruments)
    cat("\nResiduals:\n")
    print(summary(residuals(object)))
    cat("\n")
    df <- object$n - object$p
    std.errors <- sqrt(diag(object$V))
    b <- object$coefficients
    t <- b/std.errors
    p <- 2*(1 - pt(abs(t), df))
    table <- cbind(b, std.errors, t, p)
    rownames(table) <- names(b)
    colnames(table) <- c("Estimate","Std. Error","t value","Pr(>|t|)")
    print(table)
    cat(paste("\nResidual standard error:", round(object$s, digits),
        "on", df, "degrees of freedom\n\n"))
    }
    
residuals.tsls <- function(object, ...){
    res <- object$residuals
    if (is.null(object$na.action)) 
        res
    else naresid(object$na.action, res)
    }

coefficients.tsls <- function(object, ...){
    object$coefficients
    }
    
fitted.tsls <- function(object, ...){
    yhat <- as.vector(object$model.matrix %*% object$coefficients)
    if (is.null(object$na.action)) 
        yhat
    else napredict(object$na.action, yhat)
    }
    
anova.tsls <- function(object, model.2, s2, dfe, ...){
    if(class(model.2) != "tsls") stop('requires two models of class tsls')
    s2.1 <- object$s^2
    n.1 <- object$n 
    p.1 <- object$p
    dfe.1 <- n.1 - p.1
    s2.2 <- model.2$s^2
    n.2 <- model.2$n
    p.2 <- model.2$p
    dfe.2 <- n.2 - p.2
    SS.1 <- s2.1 * dfe.1
    SS.2 <- s2.2 * dfe.2
    SS <- abs(SS.1 - SS.2)
    Df <- abs(dfe.2 - dfe.1)
    if (missing(s2)){
        s2 <- if (dfe.1 > dfe.2) s2.1 else s2.2
        f <- (SS/Df) / s2
        RSS <- c(SS.1, SS.2)
        Res.Df <- c(dfe.1, dfe.2)
        SS <- c(NA, SS)
        P <- c(NA, 1 - pf(f, Df, min(dfe.1, dfe.2)))
        Df <- c(NA, Df)
        f <- c(NA, f)
        rows <- c("Model 1", "Model 2")
        }
    else{
        f <- (SS/Df) / s2
        RSS <- c(SS.1, SS.2, s2*dfe)
        Res.Df <- c(dfe.1, dfe.2, dfe)
        SS <- c(NA, SS, NA)
        P <- c(NA, 1 - pf(f, Df, min(dfe.1, dfe.2)), NA)
        Df <- c(NA, Df, NA)
        f <- c(NA, f, NA)
        rows <- c("Model 1", "Model 2", "Error")
        }
    table <- data.frame(Res.Df, RSS, Df, SS, f, P)
    head.1 <- paste("Model 1: ",format(object$formula), "  Instruments:", 
        format(object$instruments))
    head.2 <- paste("Model 2: ",format(model.2$formula), "  Instruments:", 
        format(model.2$instruments))
    names(table) <- c("Res.Df", "RSS", "Df", "Sum of Sq", "F", "Pr(>F)")
    row.names(table) <- rows
    structure(table, heading = c("Analysis of Variance", "", head.1, head.2, ""), 
        class = c("anova", "data.frame"))
    }
