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