1
$\begingroup$

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[[1]]
mu_2_hat<-mixture_model$mu[[2]]
mu_3_hat<-mixture_model$mu[[3]]
mu_4_hat<-mixture_model$mu[[4]]
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[7])
  }
  rm(m)

  #MC samples for eta
  for (m in 1:s){
    mc_samples_eta[,,m]<-(theta_k[1]*(mc_samples_xi[,1,m]))+(theta_k[2]*(mc_samples_xi[,2,m]))+(theta_k[3]*(mc_samples_xi[,3,m]))+(theta_k[4]*(mc_samples_xi[,1,m]*mc_samples_xi[,2,m]))+(theta_k[5]*(mc_samples_xi[,1,m]*mc_samples_xi[,3,m]))+(theta_k[6]*(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[7])

    }
  }
  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)
 }
$\endgroup$
3
  • $\begingroup$ Can you provide a link or a better reference to the papers you cite? $\endgroup$ – Jeremy Miles Aug 19 '17 at 18:51
  • 1
    $\begingroup$ I have provided a link to the full article. Thanks $\endgroup$ – Linniker Aug 21 '17 at 15:14
  • $\begingroup$ Any feedback would be highly appreciated $\endgroup$ – Linniker Aug 21 '17 at 16:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.