This R program is for ridge GLS estimation as in the paper \u201cOptimizing Ridge Generalized Least Squares for Structural Equation Modeling.”<\/p>\n
Authors: Miao Yang and Ke-Hai Yuan<\/p>\n
<\/p>\n
#—————————————————————————–
\n# R functions for ridge GLS estimation (don’t change codes in the functions)
\n#—————————————————————————–
\n# Computing the value for the ridge tuning parameter
\na_pred<-function(N,p,m,q,df,sk,kt){
\nm_sqrt <- sqrt(m)
\nm_2 <- m^2
\nm_log_N_2 <- log(m)\/sqrt(N)
\nm_N <- m\/N
\nq_log <- log(q)
\nq_2_N <- q^2\/(N)
\ndf_log <- log(df)
\ndf_2 <- df^2
\ndf_log_logN <- log(df)\/log(N)
\ndf_log_N <- log(df)\/(N)
\ndf_N <- df\/N
\ndf_3_2_N <- (df)^(3\/2.0)\/(N)
\nsk_3_2 <- (sk)^(3\/2.0)
\nsk_N <- sk\/N
\nkt_log <- log(kt)
\nkt_log_logN <- log(kt)\/log(N)
\nkt_sqrt_logN <- sqrt(kt)\/log(N)
\nkt_log_N_2 <- log(kt)\/sqrt(N)
\nkt_log_N <- log(kt)\/(N)
\nkt_sqrt_N <- sqrt(kt)\/(N)<\/p>\n
cb= -25.5726+0.0006*N-1.9260*m_sqrt-0.0049*m_2+24.2975*m_log_N_2-28.7679*m_N-
\n0.0331*q_log+0.0044*q_2_N+2.2866*df_log-(9.490000e-08)*df_2+8.3315*df_log_logN+
\n11.8695*df_log_N-0.0045*df_N-0.0010*df_3_2_N+0.0070*sk_3_2-732.8503*sk_N+
\n5.8133*kt_log+73.4353*kt_log_logN+0.0150*kt_sqrt_logN+192.3647*kt_log_N_2+
\n378.9057*kt_log_N+760.7915*kt_sqrt_N
\npi <- exp(cb)\/(1+exp(cb))<\/p>\n
if (pi>0.5) {
\na <- 0.05
\n}else if (pi<0.5){
\ny <- 2.7380-3.9151*log(df)\/log(N)-3.1905*log(kt)+6.9506*kt^2\/sqrt(N)
\na <- exp(y)\/(1+exp(y))
\n}<\/p>\n
return(a)
\n}<\/p>\n
# Computing relative multivariate skewness and kurtosis
\nMardia.tran1 <- function(x){
\nn <- nrow(x)
\np <- ncol(x)
\nx.c <- x-matrix(apply(x,2,mean),nrow=n,ncol=p,byrow=T)<\/p>\n
S <- cov(x)*(n-1)\/n
\nS_inv <- solve(S)
\nD <- as.matrix(x.c)%*%S_inv%*%t(x.c)<\/p>\n
b1 <- sum(D^3)
\nb2 <- sum(diag(D^2))\/n #Mardia’s kurtosis<\/p>\n
df <- p*(p+1)*(p+2)
\nsk <- b1\/(n*df)
\nkt <- b2\/(p*(p+2))<\/p>\n
return(list(‘skew’=sk,’kurt’=kt))
\n}<\/p>\n
# The following is to produce the matrix Gamma (Sy)
\nSy <- function(n,p,x){
\nps <- p*(p+1)\/2
\nmean_xt <- apply(x,2,mean)
\ny <- matrix(0,n,ps)
\nx_c <- x-matrix(1,n,1)%*%mean_xt
\nx_ct <- t(x_c)
\nScov <- x_ct %*% x_c \/n
\nvscov <- matrixcalc::vech(Scov)
\nfor (i in 1:n){
\nl <- 0
\nfor (j in 1:p){
\nfor (k in j:p){
\nl <- l+1
\ny[i,l] <- x_c[i,k]*x_c[i,j]
\n}
\n}
\n}
\nvscov_t <- t(vscov)
\ny <- y-matrix(1,n,1)%*%vscov_t # dim: n by p(p+1)\/2
\nS_y <- t(y)%*%y\/n
\nout <- list(‘S_y’=S_y, ‘vscov’=vscov, ‘Scov’=Scov)
\nreturn(out)
\n}<\/p>\n
# Evaluate structured sigma and (dsigma\/dtheta) for CFA models (needs modification if other models are used)
\nsigdsig <- function(p,m,theta0){
\nps <- p*(p+1)\/2
\np_h <- p\/m
\np2 <- p*2
\nlamb <- matrix(0,p,m)
\nfor (j in 1:m){
\np_hj <- (j-1)*p_h+1
\np_hj1 <- j*p_h
\nlamb[p_hj:p_hj1,j] <- theta0[p_hj:p_hj1]
\n}<\/p>\n
psi_vec <- theta0[(p+1):p2]
\nphi <- matrix(1,m,m)
\nk <- 0
\nfor (i in 1:(m-1)){
\nfor (j in (i+1):m){
\nk <- k+1
\nphi[j,i] <- theta0[p2+k]
\nphi[i,j] <- theta0[p2+k]
\n}
\n}
\nSigma0 <- lamb%*%phi%*%t(lamb) + diag(psi_vec)<\/p>\n
# The following is to evaluate (dSigma0\/dtheta)
\n# DL is the derivative with Lambda
\nDL <- matrix(0,ps,p)
\nlambphi <- lamb%*%phi
\nlambphi_t <- t(lambphi)
\nfor (j in 1:m){
\np_hj <- (j-1)*p_h+1
\np_hj1 <- j*p_h
\nfor (k in p_hj:p_hj1){
\ntt <- matrix(0,p,m)
\ntt[k,j] <- 1
\nttt <- tt%*%lambphi_t+lambphi%*%t(tt)
\nDL[,k] <- matrixcalc::vech(ttt)
\n}
\n}
\n#Dps is the derivative with Psi
\nDps <- matrix(0,ps,p)
\nfor (j in 1:p){
\ntt <- matrix(0,p,p)
\ntt[j,j] <- 1
\nDps[,j] <- matrixcalc::vech(tt)
\n}
\n#Dphi is the derivative with phi
\nms <- m*(m-1)\/2
\nDphi <- matrix(0,ps,ms)
\nk <- 0
\nfor (i in 1:(m-1)){
\nfor (j in (i+1):m){
\nk <- k+1
\ntt <- matrix(0,m,m)
\ntt[j,i] <- 1 ; tt[i,j] <- 1
\nttt <- lamb%*%tt%*%t(lamb)
\nDphi[,k] <- matrixcalc::vech(ttt)
\n}
\n}
\nvdsig <- cbind(DL,Dps,Dphi)
\nout <- list(‘Sigma0’=Sigma0, ‘vdsig’=vdsig)
\nreturn(out)
\n}<\/p>\n
# minimize Ridge GLS\/ADF using Fisher-scoring
\nminGLS <- function(n,p,m,theta0){
\nrequire(lavaan)
\nrequire(MASS)
\nrequire(matrixcalc)<\/p>\n
run.Sy <- Sy(n,p,x=as.matrix(variables))
\nS_y <- run.Sy$S_y
\nvscov <- run.Sy$vscov
\nweight <- a*S_y +(1-a)*diag(ps)<\/p>\n
ep <- 0.0001
\nps <- p*(p+1)\/2
\nq <- length(theta0)
\ndf <- ps-q
\ndiverg <- 1
\nestimate_GLS <- NA
\ntest_GLS <- NA
\ndelta <- 0.1
\nfor(i in 1:300){
\nsig <- sigdsig(p,m, theta0) # gives a list of 2: sig$Sigma0; sig$vdsig
\nvsig0 <- matrixcalc::vech(sig$Sigma0)
\ndswe <- t(sig$vdsig) %*% weight
\ndwd <- dswe %*% sig$vdsig
\nstdi <- try(solve(dwd),silent=T)
\nif(is.character(stdi)){
\ndiverg <- 1 # not converge
\nreturn(list(test_GLS=NA, estimate_GLS=NA, diverg=diverg))
\n}
\neresdu <- vscov-vsig0
\ndtheta <- stdi %*% dswe %*% eresdu
\ntheta0 <- theta0 + dtheta
\ndelta <- max(abs(dtheta))
\nif(delta<=ep) break;
\n}
\nif(i<300){
\ndiverg <- 0 #converge within 300 iteration
\n# Test start here;
\nf_rGLS <- t(eresdu) %*% weight %*% eresdu
\nT_GLS <- (n-1)*f_rGLS
\nVar <- stdi %*%dswe %*%S_y %*%t(dswe) %*%stdi
\nSE <- sqrt(diag(Var)\/(n-1))
\n# print(“theta_rGLS–SE_rGLS=”)
\nVmat <- weight-t(dswe)%*%stdi%*%dswe
\nVS_y <- Vmat %*% S_y
\ntVS_y <- sum(diag(VS_y))<\/p>\n
VS_y2 <- VS_y %*% VS_y
\ntVS_y2 <- sum(diag(VS_y2))
\ncst <- tVS_y\/df
\nT_GLS1 <- T_GLS\/cst
\n# pra1 <- 1-pchisq(T_GLS1,df)
\ncst_1 <- tVS_y2\/tVS_y
\ncst_2 <- (tVS_y*tVS_y)\/tVS_y2
\nT_GLS2 <- T_GLS\/cst_1
\ntest_GLS <- cbind(T_GLS,T_GLS1,T_GLS2)
\nestimate_GLS <- cbind(theta0,SE)<\/p>\n
}
\ncolnames(test_GLS)<-c(“T”,”T_rescaled”,”T_adjusted”)
\nreturn(list(test_GLS=test_GLS, estimate_GLS=estimate_GLS, diverg=diverg))
\n}<\/p>\n
#——————————————————–
\n# Codes that need user’s modification
\n#——————————————————–<\/p>\n
data(HolzingerSwineford1939) # a dataset for illustration
\nvariables <- HolzingerSwineford1939[,7:15]
\nn <- nrow(variables) # sample size
\np <- ncol(variables) # number of variables
\nm <- 3 # number of latent factors
\nq <- 2*p+m*(m-1)\/2 # number of parameters
\nps <- p*(p+1)\/2
\ndf <- ps-q # degrees of freedom
\nsk <- Mardia.tran1(variables)$skew #skewness
\nkt <- Mardia.tran1(variables)$kurt #kurtorsis
\na <- a_pred(N=n,p,m,q,df,sk,kt) #predicted value for the optimal ridge tuning paramter<\/p>\n
theta0 <- rep(1,q) #starting values<\/p>\n
minGLS(n,p,m,theta0) # output includes a list of test statistics,
\n# parameter estimates by ridge GLS and the corresponding SEs,
\n# and whether converged or not (diverge=0 means converged)<\/p>\n
<\/p>\n