search for: hpdinterval

Hi all R mailing listers: I am using the coda package. I tried to view the source of HPDinterval code by typing fix(HPDinterval), it dispalys as follows: function (obj, prob = 0.95, ...) UseMethod("HPDinterval") Then I search the answers about this case (see below), it still failed. Thank you in advance! David > getAnywhere('HPDinterval')2 differing objects matching ‘...

...mlm <- lmer(Score ~ division + (1|school), data=Age6m) (To those of you to whom this model looks familiar, thanks for your patience with this & my other questions.) Anyway, I was trying this to look at the significance of my fixed effects: A6post <- mcmcsamp(A6mlm, 50000) library(coda) HPDinterval(A6post) ..but I got this message: "no applicable method for 'HPDinterval' applied to an object of class "merMCMC"" Should I be coercing A6post to another type, or am I missing other steps altogether? Thanks :) Doug Adams ----- Doug Adams MStat Student University of...

Hi, Can anyone tell me why I am not getting the correct intervals for fixed effect terms for the following generalized linear mixed model from HPDinterval: > sessionInfo() R version 2.4.1 (2006-12-18) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] "stats" "graphi...

...advance! > data(ratdrink, package = 'faraway') > rd.er <- lmer(wt ~ weeks*treat + (1 | subject), data = ratdrink) > rd.mc <- mcmcsamp(rd.er, 10000) > library(coda) Loading required package: lattice > summary(rd.mc) Length Class Mode 1 merMCMC S4 > HPDinterval(rd.mc) Error in UseMethod("HPDinterval") : no applicable method for "HPDinterval" > str(rd.mc) Formal class 'merMCMC' [package "lme4"] with 9 slots ..@ Gp : int [1:2] 0 27 ..@ ST : num [1, 1:10000] 1.179 0.878 0.864 0.760 0.614 ... ..@ cal...

...s) samp <- rnorm(n=10000) mcmcpvalue <- function(samp) {std <- backsolve(chol(var(samp)), cbind(0,t(samp)) - colMeans(samp), transpose = TRUE) sqdist <- colSums(std*std) sum(sqdist[-1] > sqdist[1]/nrow(samp) } markov1 <- mcmcsamp(lnmass, 10000) HPDinterval(markov1) mcmcpvalue(as.matrix(markov1[,1])) mcmcpvalue(as.matrix(markov1[,2])) mcmcpvalue(as.matrix(markov1[,3])) mcmcpvalue(as.matrix(markov1[,4])) mcmcpvalue(as.matrix(markov1[,5])) mcmcpvalue(as.matrix(markov1[,6])) The first time I tried it HPDinterval generated CIs for each treatment level,...

...with Coda, I used the command, coda.options(combine.plots=TRUE, combine.stats=TRUE) to ask for results that combine the four separate chains. Sure enough, if I enter "summary(coda_odp_gini_only)", I am given a single set of output combining the four chains. However, if I enter "HPDinterval(coda_odp_gini_only)" I receive 4 sets of HPD figures, one for each chain. Is it possible to combine the four chains together to receive a single set of HPD estimates? In a similar vein, is it possible to use the Coda object to estimate the proportion of a given parameters distribution which...

...Fixed Effects: (Intr) weeks trtthr trtthy wks:trtthr weeks -0.448 treatthircl -0.707 0.317 treatthyrxn -0.642 0.288 0.454 wks:trtthrc 0.317 -0.707 -0.448 -0.203 wks:trtthyr 0.288 -0.642 -0.203 -0.448 0.454 > rd.mc <- mcmcsamp(rd.er, 50000) > library(coda) > HPDinterval(rd.mc) lower upper (Intercept) 46.420404 59.406398 weeks 25.070131 27.930363 treatthiouracil -4.420942 14.009291 treatthyroxine -10.758369 9.435761 weeks:treatthiouracil -11.404620 -7.337025 weeks:treatthyroxine -1....

...pulation genetics) could be estimated using the output of mcmcsamp (the series on mcmc sample estimates of variance components). What I have started to do is to use the matrix output that included the log(variances), exponentiate, calculate the relevant ratio, and apply either quantile or or HPDinterval to get confidence intervals. This seems too simple but I can't think of what is wrong with it. All thoughts appreciated. -Hank Dr. Hank Stevens, Assistant Professor 338 Pearson Hall Botany Department Miami University Oxford, OH 45056 Office: (513) 529-4206 Lab: (513) 529-4262 FAX: (513) 5...

Dear I am trying to find a 90% credible interval. I am using the following code. fit<-bugs( model.file=BUGScode, data=data, inits = list(geninits1,geninits2), parameters.to.save=keepers, n.chains=nchains, n.iter=runs, n.burnin=burn, n.thin=nthin, DIC= TRUE, bugs.directory="C:/Program Files/WINBUGS.14", \ ) But this is only giving 95%

Hi! I am trying to calculate HPD for the coeficients of regression models fitted with lm or lmrob in R, pretty much in the same way that can be accomplished by the association of mcmcsamp and HPDinterval functions for multilevel models fitted with lmer. Can anyone point me in the right direction on which packages/how to implement this? Thanks for your time! R. [[alternative HTML version deleted]]

...0.737 T4 -0.577 0.570 0.430 > The p-values from mcmc are: > mcmcpvalue<-function(samp) { std<-backsolve(chol(var(samp)), cbind(0,t(samp))-colMeans(samp), transpose=TRUE) sqdist<-colSums(std*std) sum(sqdist[-1]>sqdist[1]/nrow(samp) } fitSI<-mcmcsamp(fit,50000) library(coda) HPDinterval(fitSI) lower upper Intercept -4.0778905 -3.1366836 Treatment2 3.4455972 4.3196598 Treatment 3 0.399302 1.287747 Treatment 4 -1.7898933 -0.2980325 log(Treatment*Site.(in)) -22.2198233 -19.7342530 log(Site.(In)) -28.7857601 -23.0952939 mcmcpvalue(as.Matrix(fitSI[,1]))...

library(lme4) library(coda) library(languageR) fit = lmer(Reaction~Days + (1|Subject) + (0+Days|Subject), data=sleepstudy) pvals.fnc(fit)$random # compare with... samp = mcmcsamp(fit, n=10000, trans=FALSE) HPDinterval(samp) densityplot(samp, plot=F) # 'pvals.fnc' reports sigma instead of sigma^2, but it looks like the # Sbjc.(In) and Sbjc.Days are also sqrt compared with the corresponding # results from HPDinterval and densityplot. Is that intended? best wishes, Mike Allerhand

...?: det Loading required package: lme4 Attaching package: ?lme4? The following object(s) are masked from ?package:stats?: AIC, BIC Loading required package: R2WinBUGS Loading required package: coda Attaching package: ?coda? The following object(s) are masked from ?package:lme4?: HPDinterval Loading required package: abind Loading required package: foreign Attaching package: ?arm? The following object(s) are masked from ?package:coda?: traceplot ---------------------------------------------------------------------- - output with complete data ----------------------------------...

...0.1091 0.0491 2.2228 loghab 0.0875 0.0732 1.1954 landscape_diversity 1.0234 0.4850 2.1099 Correlation of Fixed Effects: (Intr) lgptch loghab logpatch 0.091 loghab -0.637 -0.121 lndscp_dvrs -0.483 -0.098 -0.348 markov1=mcmcsamp(m1,5000) HPDinterval(markov1) mcmcpvalue(as.matrix(markov1)[,1]) Error in colMeans(samp) : 'x' must be an array of at least two dimensions

...0.737 T4 -0.577 0.570 0.430 > The p-values from mcmc are: > mcmcpvalue<-function(samp) { std<-backsolve(chol(var(samp)), cbind(0,t(samp))-colMeans(samp), transpose=TRUE) sqdist<-colSums(std*std) sum(sqdist[-1]>sqdist[1]/nrow(samp) } fitSI<-mcmcsamp(fit,50000) library(coda) HPDinterval(fitSI) lower upper Intercept -4.0778905 -3.1366836 Treatment2 3.4455972 4.3196598 Treatment 3 0.399302 1.287747 Treatment 4 -1.7898933 -0.2980325 log(Treatment*Site.(in)) -22.2198233 -19.7342530 log(Site.(In)) -28.7857601 -23.0952939 mcmcpvalue(as.Matrix(fitSI[,1])) e...

Dear all, I am running a mixed model using lmer. In order to obtain CI of individual coefficients I use mcmcsamp. However, I need advice which values that are most appropriate to present in result section of a paper. I have not used mixed models and lmer so much before so my question is probably very naive. However, to avoid to much problems with journal editors and referees addicted to

...BIC logLik Chisq Chi Df Pr(>Chisq) m2c 4 89.036 96.173 -40.518 m2 5 84.834 93.755 -37.417 6.2023 1 0.01276 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 The p-values from mcmc are: ## markov1=mcmcsamp(m2,5000) HPDinterval(markov1) lower upper (Intercept) -1.394287660 0.6023229 logpatch 0.031154910 0.1906861 loghab 0.002961281 0.2165285 landscape_diversity 0.245623183 1.6442544 log(site.(In)) -41.156007604 -1.6993996 attr(,"Probabi...

...CE and REMOVE fns (Tina Robles) > 48. Re: (Newbie)Basic Basic global vs. local > variables (projection83) > 49. Re: Reading user input (projection83) > 50. Newbie: Simple loops: complex troubles > (projection83) > 51. gnomeGUI under KDE (Luis Lopez Oliveros) > 52. Re: HPDinterval problem (Seyed Reza > Jafarzadeh) > 53. Re: converting a list to a data.frame (Patrick > Connolly) > 54. time zone problems (Marc Fischer) > 55. Re: Reading user input (Petr Klasterecky) > 56. Re: Newbie: Simple loops: complex troubles > (Petr Klasterecky) > 57. Re...

...lues? > > Below is the corresponding R code with some output so that you can see: > > ## > fit<-lmer(End~Treatment+offset(log(Area)+(1|Site/Treatment), family=poisson) > Summary Intercept > The p-values from mcmc are: > > ## > markov1=mcmcsamp(m2,5000) > > HPDinterval(markov1) > lower upper > (Intercept) -1.394287660 0.6023229 > logpatch 0.031154910 0.1906861 > loghab 0.002961281 0.2165285 > landscape_diversity 0.245623183 1.6442544 > log(site.(In)) -41.156007604 -1.6993996 > attr(,"Probability") > [1] 0.95 > > ## > &gt...

When I do lmer models I only get Estimate, Standard Error and t value in the output for the fixed effects. Is there a way I get degrees of freedom and p values as well? I'm a very new to R, so sorry if this a stupid question. Thank you - Mike Mike Ford Centre for Speech and Language Department of Experimental Psychology Downing Street Cambridge CB2 3EB Tel: +44 (0) 1223 766559 Fax: +44