Hi,
I have a dataset with the binary outcome Y(0,1) and 4 covariates (X1,X@,X#,X$).
I am trying to use MCMClogit to model logistic regression using MCMC. I am
getting an error where it doesnt identify the covariates ,although its reading
in correctly. The dataset is a sample of actual dataset. Below is my code:
> #######################
>
>
> #retreive data
> # considering four covariates
> d.df=as.data.frame(read.table("c:/tina/phd/thesis/data/modified_data1.1.txt",header=T,sep=","))
> y=d.df[,ncol(d.df)]
> x=d.df[,1:4]
> c.df=cbind(y,x)
> #x=cbind(1,x)
> p <- ncol(c.df)
>
> # marginal log-prior of beta[]
> logpriorfun <- function(beta, mu, gshape, grate)
+ {
+ logprior = -p*log(2) + log(gamma(p+gshape)) - log(gamma(gshape))
+ + gshape*log(grate) - (p+gshape)* log(grate+sum(abs(beta)))
+ return(logprior)
+ }
> require(MCMCpack)
Loading required package: MCMCpack
Loading required package: coda
Loading required package: lattice
Loading required package: MASS
##
## Markov Chain Monte Carlo Package (MCMCpack)
## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
##
## Support provided by the U.S. National Science Foundation
## (Grants SES-0350646 and SES-0350613)
##
[1] TRUE
Warning message:
package 'MASS' was built under R version 2.4.1
> a0 = 0.5
> b0 = 1
> mu0 = 0
> beta.init=list(c(0, rep(0.1,4)), c(0, rep(-0.1,4)), c(0, rep(0, 4)))
> burnin.cycles = 1000
> mcmc.cycles = 25000
> # three chains
> post.list <- lapply(beta.init, function(vec)
+ {
+ posterior <- MCMClogit(y~x1+x2+x3+x4, data=c.df, burnin=burnin.cycles,
mcmc=mcmc.cycles,
+ thin=5, tune=0.5, beta.start=vec, user.prior.density=logpriorfun, logfun=TRUE,
+ mu=mu0, gshape=a0, grate=b0)
+ return(posterior)
+ })
Error in eval(expr, envir, enclos) : object "x1" not found
>
Any suggestions will be greatly appreciated.
Thanks,
Anamika
---------------------------------
We won't tell. Get more on shows you hate to love
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
