REPLY TO ALL FOR THE R-HELP LIST!!!
I apologize for the bluntness but you must realize that it is critical if
you desire to get help on the mailing list beyond a single person your
question must actually get to the mailing list. My expertise only goes so
far and there is an infinitely larger community that could likely help you
more quickly and often in a way that can increase performance. I am always
happy to help when time allows but please make sure the r-help@r-project.org
is a recipient as well.
Regarding your code, after glancing at it again it is clear that you aren't
actually submitting a covariance matrix to mvrnorm. The only place you
specify phi1 is initially with NA's and then within the list init1. If you
intend to use the list init1 you specify it within the function call:
xi1<-mvrnorm(1,c(0,0,0,0),init1$phi1, tol = 1e-6, empirical = FALSE,
EISPACK = FALSE)
Once again, please make sure you are not just sending these messages solely
to me.
Regards,
Charles
On Thu, Jun 5, 2014 at 9:04 AM, thanoon younis <thanoon.youni...@gmail.com>
wrote:
> many thanks to you Dr. Charles
>
> Really i have a problem with simulation data in xi and now i have this
> erro r "Error in mvrnorm(1, c(0, 0, 0), phi1, tol = 1e-06, empirical =
> FALSE, : incompatible arguments"
>
> Regards
>
>
>
>
> On 5 June 2014 16:45, Charles Determan Jr <deter...@umn.edu> wrote:
>
>> 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.youni...@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
>>
>
>
--
Dr. Charles Determan, PhD
Integrated Biosciences
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
