[R] MCMC logit
Petr Klasterecky
klaster at karlin.mff.cuni.cz
Mon Mar 12 08:09:48 CET 2007
Please, DO read your error messages - it says
"Error in tune %*% V : non-conformable arguments" - most probably you
are trying to matrix-multiply (%*%) a scalar (tune) and a vector or
matrix (V)...
This has nothing to do with the previous error message "Error in
eval(expr, envir, enclos) : object "x1" not found" - this one was fixed
by introducing separate covariates x1-x4. Btw, aparently you DO NOT know
that you dataframe has names exactly as you want them - R is
case-sensitive and thus x1 is NOT the same as X1... However, you are
supposed to read R-intro, R-FAQ etc before posting to the list, so you
already know all this...
Petr
Anamika Chaudhuri napsal(a):
> Hi,
>
> I know the data frame c.df has the variables named exactly the way I want it to be. I tried reading each covariate seperately but still it gives me an error. I tried fidning out the error but it doesnt seem very evident. Here is the error message
>
> > #######################
>>
>> #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]
>> # read each column of x seperately
>> x1=d.df[,1]
>> x2=d.df[,2]
>> x3=d.df[,3]
>> x4=d.df[,4]
>> c.df=cbind(y,x)
>> print(c.df)
> y X1 X2 X3 X4
> 1 1 2 0 0 1
> 2 1 10 0 0 1
> 3 0 4 0 0 1
> 4 1 0 0 1 1
> 5 1 0 0 1 1
> 6 1 7 0 0 1
> 7 1 1 0 0 1
> 8 0 0 0 0 1
> 9 1 0 0 0 1
> 10 1 5 0 0 1
>> p <- ncol(c.df)
>>
>> # setting error handler to true
>> options(show.error.mesages = TRUE)
>>
>> # 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 tune %*% V : non-conformable arguments
>> # for tracing error mesages
>> geterrmessage()
> [1] "Error in tune %*% V : non-conformable arguments\n"
>> traceback()
> 3: 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)
> 2: FUN(X[[1]], ...)
> 1: 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)
> })
>
>
> Thanks for your help,
> Anamika
> Ravi Varadhan <rvaradhan at jhmi.edu> wrote:
> As the error message clearly indicates, the function MCMClogit is unable to
> find the variable x1 (possibly x2,x3, and x4 also) in the data frame c.df.
> Check the names of the variables in that data frame and make sure that the
> names correspond to the formula specification.
>
> Hope this helps,
> Ravi.
>
> ----------------------------------------------------------------------------
> -------
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: rvaradhan at jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----------------------------------------------------------------------------
> --------
>
>
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Anamika Chaudhuri
> Sent: Friday, March 09, 2007 3:27 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] MCMC logit
>
> 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]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>
>
>
> ---------------------------------
> Bored stiff? Loosen up...
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>
--
Petr Klasterecky
Dept. of Probability and Statistics
Charles University in Prague
Czech Republic
More information about the R-help
mailing list