1cluster.post.hook <- function(x,...) {
2    if (class(x)[1]=="multigroupfit") {
3        if (is.null(x$cluster)) return(NULL)
4        if (any(unlist(lapply(x$cluster,is.null)))) return(NULL)
5        allclusters <- unlist(x$cluster)
6        uclust <- unique(allclusters)
7        K <- length(uclust)
8        G <- x$model$ngroup
9        S0 <- lapply(score(x,indiv=TRUE), function(x) { x[which(is.na(x))] <- 0; x })
10        S <- matrix(0,length(pars(x)),nrow=K)
11        for (i in uclust) {
12            for (j in seq_len(G)) {
13                idx <- which(x$cluster[[j]]==i)
14                if (length(idx)>0)
15                    S[i,] <- S[i,] + colSums(S0[[j]][idx,,drop=FALSE])
16            }
17        }
18        J <- crossprod(S)
19        I <- information(x,type="hessian",...)
20        iI <- Inverse(I)
21        asVar <- iI%*%J%*%iI
22        x$vcov <- asVar
23        return(x)
24    }
25
26    ## lvmfit:
27    if (!is.null(x$cluster)) {
28        uclust <- unique(x$cluster)
29        K <- length(uclust)
30        S <- score(x,indiv=TRUE) #,...)
31        I <- information(x,type="hessian") #,...)
32        iI <- Inverse(I)
33
34        S0 <- matrix(0,ncol=ncol(S),nrow=K)
35        count <- 0
36        for (i in uclust) {
37            count <- count+1
38            S0[count,] <- colSums(S[which(x$cluster==i),,drop=FALSE])
39        };
40        adj1 <- K/(K-1)
41        ## p <- ncol(S)
42        ## adj1 <- K/(K-p) ## Mancl & DeRouen, 2001
43
44        J <- adj1*crossprod(S0)
45        col3 <- sqrt(diag(iI)); ## Naive se
46        nn <- c("Estimate","Robust SE", "Naive SE", "P-value")
47        asVar <- iI%*%J%*%iI
48    } else {
49        asVar <- x$vcov
50    }
51    diag(asVar)[diag(asVar)==0] <- NA
52    mycoef <- x$opt$estimate
53    x$vcov <- asVar
54    SD <- sqrt(diag(asVar))
55    Z <- mycoef/SD
56    pval <- 2*(pnorm(abs(Z),lower.tail=FALSE))
57    if (is.null(x$cluster)) {
58        col3 <- Z
59        nn <-  c("Estimate","Std. Error", "Z-value", "P-value")
60    }
61    newcoef <- cbind(mycoef, SD, col3, pval);
62    nparall <- index(x)$npar + ifelse(x$control$meanstructure, index(x)$npar.mean,0)
63    if (!is.null(x$expar)) {
64        nparall <- nparall+length(x$expar)
65    }
66    mycoef <- matrix(NA,nrow=nparall,ncol=4)
67    mycoef[x$pp.idx,] <- newcoef
68    colnames(mycoef) <- nn
69    mynames <- c()
70    if (x$control$meanstructure) {
71        mynames <- vars(x)[index(x)$v1==1]
72    }
73    if (index(x)$npar>0) {
74        mynames <- c(mynames,paste0("p",seq_len(index(x)$npar)))
75    }
76    if (!is.null(x$expar)) {
77        mynames <- c(mynames,names(x$expar))
78    }
79
80    rownames(mycoef) <- mynames
81    x$coef <- mycoef
82    return(x)
83}
84