# Nonlinear Structural Equation Modelling- Estimation via Pseudo Maximum Likelihood in R

I am trying to implement the pseudo maximum likelihood approach described in Wall and Amemiya (2007) (see link to article at the end) in estimating the parameters of a nonlinear structural model of the form:

$\eta=\gamma_1\xi_1+ \gamma_2\xi_2+ \gamma_3\xi_3+\omega_{12}\xi_1\xi_2+\omega_{13}\xi_1\xi_3+\omega_{23}\xi_{2}\xi_{3}+\zeta$

where $\zeta\sim N(0,\Delta)$ and the measurement model is given by:

$\boldsymbol{X}=\boldsymbol{\Lambda}_{x}\boldsymbol{\xi} + \boldsymbol{\delta}_{x}$

$\boldsymbol{Y}=\boldsymbol{\Lambda}_{y}\eta + \boldsymbol{\delta}_{y}$

which is combined into one model as:

$\boldsymbol{Z}=\boldsymbol{\Lambda}\boldsymbol{f} + \boldsymbol{\epsilon}$

where $\boldsymbol{Z}=(\boldsymbol{Y}',\boldsymbol{X}')'$ and $\boldsymbol{f}=(\eta,\boldsymbol{\xi}')'$ and $\boldsymbol{\epsilon}=( \boldsymbol{\delta}_{y}', \boldsymbol{\delta}_{x}')'$.

The method starts by first estimating the parameters in the measurement model through maximum likelihood. Using these parameter estimates, Bartlett scores for $\eta$ and $\boldsymbol{\xi}$ are computed. Using these factor scores, a Gaussian Mixture model is fitted to $\boldsymbol{\xi}$.

As described in the paper attached, the parameters in the structural model are then estimated via a Monte-Carlo EM algorithm.

I have written a script in R to perform a simulation to assess parameter recovery properties of the method. However, the algorithm is not converging. I suspect that I've some problems in the code as I am still a novice in R programming. My gut feeling is that I have a mistake in the calculation of the weights but I cannot seem to find the problem. I've been going at it for weeks on end to no avail.

I am attaching a copy of the script I have written in the hope of getting some feedback or a new perspective as to where the problem with the code might be.

Many thanks!

Paper refered to in text: Wall and Amemiya (2007)- Nonlinear Structural Equation Modeling as a Statistical Method: Can be accessed on: https://pdfs.semanticscholar.org/590f/da1c7f216a110d4c11fff2d61eee2f451462.pdf

#Load Required Packages
library("lavaan", lib.loc="~/R/win-library/3.4")
library("mclust", lib.loc="~/R/win-library/3.4")
library("mixtools", lib.loc="~/R/win-library/3.4")
library("MASS", lib.loc="~/R/win-library/3.4")
library("MASS", lib.loc="C:/Program Files/R/R-3.4.1/library")
library("MixSim", lib.loc="~/R/win-library/3.4")
library("mvtnorm", lib.loc="~/R/win-library/3.4")
library("mixAK", lib.loc="~/R/win-library/3.4")
library("gdata", lib.loc="~/R/win-library/3.4")

#Data Simulation

n=1000 #sample size
s=100 #number of Monte-Carlo samples

#Generate population latent variables
mu_z<-c(rep(0,3))
sigma_z<-matrix(rep(0,9),3,3)
sigma_z[,1]<-c(1,0.5,0.6)
sigma_z[,2]<-c(0.5,1,0.4)
sigma_z[,3]<-c(0.6,0.4,1)

#Compute xi from zs (as is done in paper)
set.seed(28)
z<-matrix(0,ncol=3,nrow=n)
z<-mvrnorm(n, mu=mu_z, Sigma=sigma_z)
xi<-sqrt(exp(z))

#Error term of structual model
var_zeta<-1
zeta<-rnorm(n,0,sqrt(var_zeta))

#Setting Population parameters
gamma1<-1.5
gamma2<-1
gamma3<-2
omega12<--3
omega13<-1.25
omega23<-1.75

eta<-(gamma1*(xi[,1]))+(gamma2*(xi[,2]))+(gamma3*(xi[,3]))+(omega12*(xi[,1]*xi[,2]))+(omega13*(xi[,1]*xi[,3]))+(omega23*(xi[,2]*xi[,3]))+(zeta)

f<-cbind(eta,xi)

#Computing observable/indicator variables

#Variance-covariance matrix of error term in measurement model
cov_delta_epsilon<-diag(0.25,9)

#Measurement error
set.seed(281)
error_measurement_model<-mvrnorm(n,mu=c(rep(0,9)),Sigma=cov_delta_epsilon,empirical=FALSE)

#Matrix of observable variables
observable<-matrix(0,ncol=9,nrow=n)
observable[,1]<-f[,1]+error_measurement_model[,1]
observable[,2]<-f[,1]+error_measurement_model[,2]
observable[,3]<-f[,1]+error_measurement_model[,3]
observable[,4]<-f[,2]+error_measurement_model[,4]
observable[,5]<-f[,2]+error_measurement_model[,5]
observable[,6]<-f[,3]+error_measurement_model[,6]
observable[,7]<-f[,3]+error_measurement_model[,7]
observable[,8]<-f[,4]+error_measurement_model[,8]
observable[,9]<-f[,4]+error_measurement_model[,9]

as.data.frame(observable)
colnames(observable) <- c("y1","y2","y3","X1","X2","X3","X4","X5","X6")

#Estimation of Nuisance Parameters of measurment model using MLE
model<-'eta1=~y1+y2+y3
xi_1=~X1+X2
xi_2=~X3+X4
xi_3=~X5+X6'

cfa_model<-cfa(model,data=observable,estimator="ML")

#Extraction of nuisance parameters (measurement model parameters)

#matrix of estimated variances of error terms
psi_hat<-lavInspect(cfa_model,"est")$theta #matrix of estimated loadings lambda_hat<-lavInspect(cfa_model,"est")$lambda

#Computation of Bartlett factor scores
f_hat<-matrix(0,nrow=n,ncol=ncol(f))

#The following gives the matrix of estimated scores
#where the ith row is the score of the ith individual using equation 19.
for (i in 1:n){
f_hat[i,]<-solve(t(lambda_hat)%*%solve(psi_hat)%*%lambda_hat)%*%t(lambda_hat)%*%solve(psi_hat)%*%observable[i,]
}
f_hat<-as.data.frame(f_hat)

#Extract xi scores
xi_scores<-as.matrix(f_hat[,2:4])

#Computation of sigma_r_xi_hat (equation 21)
sigma_r_hat<-solve(t(lambda_hat)%*%solve(psi_hat)%*%lambda_hat)

#Extracting elements from sigma_r_hat which correspond to xi
sigma_r_xi_hat<-cbind(sigma_r_hat[2:4,2],sigma_r_hat[2:4,3],sigma_r_hat[2:4,4])

#Estimation of mixture parameters of xi_hat

#First check what type of mixture fits the xi scores best
mixture<-Mclust(xi_scores)

sigma_initial<-list(diag(1,3),diag(1,3),diag(1,3),diag(1,3))
mu_initial<-list(c(rep(0,3)),c(rep(0,3)),c(rep(0,3)),c(rep(0,3)))

#Fit mixture distribution to obtain pi_hat, mu_hat and sigma_hat_xi_hat
mixture_model<-mvnormalmixEM(xi_scores,k=4,mu=mu_initial,sigma=sigma_initial,arbmean=TRUE,arbvar=FALSE)
pi_hat<-mixture_model$lambda mu_1_hat<-mixture_model$mu[]
mu_2_hat<-mixture_model$mu[] mu_3_hat<-mixture_model$mu[]
mu_4_hat<-mixture_model$mu[] mu_hat_matrix<-rbind(mu_1_hat,mu_2_hat,mu_3_hat,mu_4_hat) sigma_hat_xi_hat<-mixture_model$sigma

#Obtain sigma_hat_xi through subtraction
sigma_hat_xi<-sigma_hat_xi_hat-sigma_r_xi_hat

#Setting initial parameter estimates for MCEM

#Initial parameter estimates to compute Monte-Carlo etas
set.seed(11)
gamma1_initial<-runif(1)
set.seed(12)
gamma2_initial<-runif(1)
set.seed(13)
gamma3_initial<-runif(1)
set.seed(14)
omega12_initial<-runif(1)
set.seed(15)
omega13_initial<-runif(1)
set.seed(16)
omega23_initial<-runif(1)

#Initial estimate of variance of zeta
set.seed(211)
est_var_zeta<-runif(1)

#Setting initial parameter vectors at iteration k and k+1
theta_k<-c(gamma1_initial,gamma2_initial,gamma3_initial,omega12_initial,omega13_initial,omega23_initial,sqrt(est_var_zeta))
theta_k_plus_1<-rep(0,length(theta_k))

diff_iter<-10

#xi MC samples

#Repeating variance-covariance matrix of xi_hat
s_mc<-array(0,dim=c(3,3,4))
for (i in 1:4){
s_mc[,,i]<-sigma_hat_xi
}

mc_samples_xi<-array(0,c(n,3,s))
for (j in 1:s){
mc_samples_xi[,,j]<-simdataset(n, Pi=pi_hat, Mu=mu_hat_matrix, S=s_mc)$X } rm(j) #Combine all xis into one matrix xi_reg<-mc_samples_xi[,,1] for (i in 2:s){ xi_reg<-rbind(xi_reg,mc_samples_xi[,,i]) } rm(i) #Define model for weighted least squares xi_one<-xi_reg[,1] xi_two<-xi_reg[,2] xi_three<-xi_reg[,3] denominator_3<-array(0,dim=c(n,s)) #Repeating variance-covariance matrix of mixture distribution as function requires a matrix for each mixture distribution sigma_hat_xi_list<-list(sigma_hat_xi,sigma_hat_xi,sigma_hat_xi,sigma_hat_xi_hat) #f(xi) from multivariate mixture of normal distributions for (i in 1:n){ for (m in 1:s){ denominator_3[i,m]<-dMVNmixture(x=mc_samples_xi[i,,m], weight=pi_hat, mean=mu_hat_matrix,Sigma=sigma_hat_xi_list, log=FALSE) } } rm(i) rm(m) ######################################################################################################################################################################################################################################################################################################################################## #MCEM Algorithm while (diff_iter>0.001){ mc_samples_zeta<-array(0,dim=c(n,1,s)) mc_samples_eta<-array(0,dim=c(n,1,s)) ##Generate Monte-Carlo samples of xis #xi monte-carlo samples were here #Generate Monte-Carlo samples of xi #MC samples for zeta for (m in 1:s){ mc_samples_zeta[,,m]<-rnorm(n,0,theta_k) } rm(m) #MC samples for eta for (m in 1:s){ mc_samples_eta[,,m]<-(theta_k*(mc_samples_xi[,1,m]))+(theta_k*(mc_samples_xi[,2,m]))+(theta_k*(mc_samples_xi[,3,m]))+(theta_k*(mc_samples_xi[,1,m]*mc_samples_xi[,2,m]))+(theta_k*(mc_samples_xi[,1,m]*mc_samples_xi[,3,m]))+(theta_k*(mc_samples_xi[,2,m]*mc_samples_xi[,3,m]))+(mc_samples_zeta[,1,m]) } rm(m) f_i_hat<-array(0,dim=c(n,4,s)) for (m in 1:s){ f_i_hat[,,m]<-cbind(mc_samples_eta[,,m],mc_samples_xi[,,m]) } #E-Step: Computation of weights #Computation of numerator of weights from multivariate normal distribution numerator<-array(0,c(n,s)) for (i in 1:n){ for (m in 1:s){ numerator[i,m]<-dmvnorm(observable[i,], mean = lambda_hat%*%f_i_hat[i,,m], sigma = psi_hat, log = FALSE) } } rm(i) rm(m) #Computation of denominator of weights mc_samples_eta_hat<-array(0,dim=c(n,1,s)) for (m in 1:s){ mc_samples_eta_hat[,,m]<-mc_samples_eta[,,m]-mc_samples_zeta[,,m] } rm(m) denominator_1<-numerator denominator_2<-array(0,dim=c(n,s)) #f(eta|xi) from univariate normal distribution for (i in 1:n){ for (m in 1:s){ denominator_2[i,m]<-dnorm(mc_samples_eta[i,,m],mean=mc_samples_eta_hat[i,,m],sd=theta_k) } } rm(i) rm(m) #Denominator 3 here #Final denominator as the product of the three terms making it up denominator<-denominator_1*denominator_2*denominator_3 denominator<-apply(denominator,1,sum)/s #Final weights W_i_m<-numerator/denominator #Perform weighted Least Squares #Combine all etas into one vector/matrix eta_reg<-mc_samples_eta[,,1] for (i in 2:s){ eta_reg<-c(eta_reg,mc_samples_eta[,,i]) } rm(i) #Combine all weights into one vector reg_weight<-W_i_m[,1] for (i in 2:s){ reg_weight<-c(reg_weight,W_i_m[,i]) } rm(i) reg_weight<-reg_weight/s reg_weight[which(is.nan(reg_weight))] = NA reg_weight[which(is.infinite(reg_weight))] = NA #Perform weighted least squares wls<-lm(eta_reg~0+xi_one+xi_two+xi_three+xi_one:xi_two+xi_one:xi_three+xi_two:xi_three,na.action=na.exclude) wls$coefficients

theta_k_plus_1<-c(wls$coefficients,summary(wls)$sigma)
print(theta_k_plus_1)

diff_iter<-as.numeric(sqrt(sum((theta_k-theta_k_plus_1)^2)))
print(diff_iter)
theta_k<-theta_k_plus_1

rm(mc_samples_zeta)
rm(mc_samples_eta)
rm(numerator)
rm(denominator)
rm(denominator_1)
rm(denominator_2)
#rm(denominator_3)
rm(W_i_m)
rm(reg_weight)
rm(eta_reg)
rm(f_i_hat)
rm(mc_samples_eta_hat)
}

• Can you provide a link or a better reference to the papers you cite? – Jeremy Miles Aug 19 '17 at 18:51
• I have provided a link to the full article. Thanks – Linniker Aug 21 '17 at 15:14
• Any feedback would be highly appreciated – Linniker Aug 21 '17 at 16:13