Hello again Thanoon,
Once again, you should send these request not to me but to the r-help
list. You are far more likely to get help from the greater R community
than just me. Furthermore, it is not entirely clear where your error is.
It is courteous to provide only the code that is run up to the error and
commenting the location. I see you have a loop and all the more important
to isolate where the error is taking place and commenting it out.
My recommendations:
1. Set both t and i = 1 and try to run your loop code manually to locate
the specific function providing the error.
2. Check your data inputs for the function. The error tells you that some
data is missing or infinite. Perhaps a slight error in your data
generation?
3. See if you can address the specific instance before running the full
loops, it may be multiple instances or just one.
I hope this provides some guidance,
Charles
On Thu, Jun 5, 2014 at 3:10 AM, thanoon younis
<thanoon.younis80@gmail.com>
wrote:
> Dear Dr. Charles
> i need your help to correct the winbugs code below to estimate parameters
> of SEM by using bayesian inference for two group. I have this
error"Error
> in eigen(Sigma, symmetric = TRUE, EISPACK = EISPACK) : infinite or missing
> values in 'x'.
>
> many thanks in advance
>
> R-CODE
>
> library(MASS) #Load the MASS package
> library(R2WinBUGS) #Load the R2WinBUGS package
> library(boa) #Load the boa package
> library(coda) #Load the coda package
>
> N1<-200;N2=100; P1<-10;P2<-10
>
> phi1<-matrix(NA,nrow=4,ncol=4) #The covariance matrix of xi1
> Ro1<-matrix(NA,nrow=4,ncol=4)
> yo1<-matrix(data=NA,nrow=200,ncol=10)
>
> phi2<-matrix(NA,nrow=4,ncol=4) #The covariance matrix of xi2
> Ro2<-matrix(NA,nrow=4,ncol=4)
> yo2<-matrix(data=NA,nrow=200,ncol=10)
>
> #Matrices save the Bayesian Estimates and Standard Errors
> Eu1<-matrix(data=NA,nrow=100,ncol=10)
> SEu<-matrix(data=NA,nrow=100,ncol=10)
> Elam1<-matrix(data=NA,nrow=100,ncol=6)
> SElam1<-matrix(data=NA,nrow=100,ncol=6)
> Egam1<-matrix(data=NA,nrow=100,ncol=4)
> SEgam1<-matrix(data=NA,nrow=100,ncol=4)
> Ephx1<-matrix(data=NA,nrow=100,ncol=4)
> SEphx1<-matrix(data=NA,nrow=100,ncol=4)
> Eb1<-numeric(100); SEb<-numeric(100)
> Esgd1<-numeric(100); SEsgd<-numeric(100)
> Eu2<-matrix(data=NA,nrow=100,ncol=10)
> SEu2<-matrix(data=NA,nrow=100,ncol=10)
> Elam2<-matrix(data=NA,nrow=100,ncol=6)
> SElam2<-matrix(data=NA,nrow=100,ncol=6)
> Egam2<-matrix(data=NA,nrow=100,ncol=3)
> SEgam2<-matrix(data=NA,nrow=100,ncol=3)
> Ephx2<-matrix(data=NA,nrow=100,ncol=3)
> SEphx2<-matrix(data=NA,nrow=100,ncol=3)
> Eb2<-numeric(100); SEb2<-numeric(100)
> Esgd2<-numeric(100); SEsgd2<-numeric(100)
>
>
> #Arrays save the HPD intervals
> mu.y1=array(NA, c(100,10,2))
> lam1=array(NA, c(100,6,2))
> gam1=array(NA, c(100,4,2))
> sgd1=array(NA, c(100,2))
> phx1=array(NA, c(100,4,2))
> mu.y2=array(NA, c(100,10,2))
> lam2=array(NA, c(100,6,2))
> gam2=array(NA, c(100,3,2))
> sgd2=array(NA, c(100,2))
> phx2=array(NA, c(100,3,2))
>
> DIC=numeric(100) #DIC values
>
> #Parameters to be estimated
> parameters<-c
>
("mu.y1","lam1","gam1","phi1","psi1","psd1","sgd1","phx1","mu.y2","lam2","gam2","phi2","psi2","psd2","sgd2","phx2")
>
> #Initial values for the MCMC in WinBUGS
> init1<-list(uby1=rep(0.0,10),lam1=rep(0.0,10),
>
>
gam1=c(1.0,1.0,1.0,1.0),psd1=1.0,phi1=matrix(data=c(1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0),ncol=4,byrow=TRUE),psi1=rep(0.0,10),
> xi1=matrix(data=rep(0.0,200),ncol=4),xi2=matrix(data=rep(0.0,200),ncol=4))
>
>
> init2<-list(mu.y1=rep(0.0,10),lam1=rep(0.0,10),
>
gam1=c(1.0,1.0,1.0,1.0),psd1=1.0,phi1=matrix(data=c(1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0),ncol=4,byrow=TRUE),psi1=rep(0.0,10),xi1=matrix(data=rep
> (0.0,200),ncol=4),xi2=matrix(data=rep(0.0,200),ncol=4))
>
>
> inits<-list(init1, init2)
>
> #Do simulation for 100 replications
> for (t in 1:100) {
> #Generate the data for the simulation study
> for (i in 1:N1) {
> #Generate xi1
> xi1<-mvrnorm(1,c(0,0,0,0),phi1, tol = 1e-6, empirical = FALSE,
> EISPACK = FALSE)
> #Generate error term is structural equation
> delta<-rnorm(1,0,sqrt(0.3))
> #Generate eta1 according to the structural equation
> eta1<-xi1[1]+xi1[2]+xi1[3]+xi1[1]*xi1[2]+delta
> #Generate error terms in measurement equations
> eps1<-rnorm(3,0,1)
>
> #Generate theta according to measurement equations
>
>
> Y1[i,1]=Xi1[i,1]+eps1[1]
> Y1[i,2]=0.8*Xi1[i,1]+eps1[2]
> Y1[i,3]=0.8*Xi1[i,1]+eps1[3]
> Y1[i,4]=0.8*Xi1[i,1]+eps1[4]
> Y1[i,5]=Xi1[i,2]+eps1[5]
> Y1[i,6]=0.8*Xi1[i,2]+eps1[6]
> Y1[i,7]=Xi1[i,3]+eps1[7]
> Y1[i,8]=Xi1[i,2]+eps1[8]
> Y1[i,9]=Eta1[i]+eps1[9]
> Y1[i,10]=0.8*Eta1[i]+eps1[10]
> }
>
> #transform theta to orinal variables
> for (j in 1:10) { if (Y1[j]>0) yo1[i,j]<-1 else
yo1[i,j]<-0 }
>
>
> #Input data set for WinBUGS
> data<-list(N1=200,N2=200,P=10,R=Ro,z=yo1)
>
> #Call WinBUGS
>
model<-bugs(data,inits,parameters,model.file="D:/Run/model.txt",
> n.chains=2,n.iter=5000,n.burnin=1,n.thin=1,DIC=True,
> bugs.directory="c:/Program Files/WinBUGS14/",
> working.directory="D:/Run/")
>
> #Save Bayesian Estimates
> Eu1[t,]<-model$mean$uby; Elam1[t,]<-model$mean$lam;
> Egam1[t,]<-model$mean$gam
> Ephx1[t,1]<-model$mean$phx1[1,1];
Ephx1[t,2]<-model$mean$phx1[1,2]
> Ephx1[t,3]<-model$mean$phx1[2,2];
> Esgd1[t]<-model$mean$sgd
>
> #Save Standard Errors
> SEu1[t,]<-model$sd$uby1; SElam1[t,]<-model$sd$lam;
> SEgam1[t,]<-model$sd$gam1
> SEphx1[t,1]<-model$sd$phx1[1,1]; SEphx1[t,2]<-model$sd$phx1[1,2]
> SEphx1[t,3]<-model$sd$phx1[2,2]; SEb1[t]<-model$sd$ubeta1
> SEsgd1[t]<-model$sd$sgd
>
> #Save HPD intervals
> for (i in 1:10) {
> temp=model$sims.array[,1,i];
> uby[t,i,]=boa.hpd(temp,0.05)
> }
> temp=model$sims.array[,1,10]; ubeta[t,]=boa.hpd(temp,0.05)
> for (i in 1:6) {
> temp=model$sims.array[,1,10+i];
> lam[t,i,]=boa.hpd(temp,0.05)
> }
> for (i in 1:3) {
> temp=model$sims.array[,1,16+i];
> gam[t,i,]=boa.hpd(temp,0.05)
> }
> temp=model$sims.array[,1,20]; sgd[t,]=boa.hpd(temp,0.05)
> temp=model$sims.array[,1,21]; phx[t,1,]=boa.hpd(temp,0.05)
> temp=model$sims.array[,1,22]; phx[t,2,]=boa.hpd(temp,0.05)
> temp=model$sims.array[,1,24]; phx[t,3,]=boa.hpd(temp,0.05)
>
> #Save DIC value
> DIC[t]=model#DIC
> } #end
>
>
>
--
Dr. Charles Determan, PhD
Integrated Biosciences
University of Minnesota
[[alternative HTML version deleted]]