[R] lmer and estimation of p-values: error with mcmcpvalue()

2007-02-12 Thread Christoph Scherber 
Dear all,

I am currently analyzing count data from a hierarchical design, and I?ve 
tried to follow the suggestions for a correct estimation of p-values as 
discusssed at R-Wiki 
(http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-testss=lme%20and%20aov).

However, I have the problem that my model only consists of parameters 
with just 1 d.f. (intercepts, slopes), so that the mcmcpvalue function 
defined below obviously produces error messages.

How can I proceed in estimating the p-values, then?

I very much acknowledge any suggestions.

Best regards
Christoph.

##
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) }

m1-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),quasipoisson)

Generalized linear mixed model fit using Laplace
Formula: number_pollinators ~ logpatch + loghab + landscape_diversity + 
 (1 | site)
Data: primula
  Family: quasipoisson(log link)
AIC   BIC logLik deviance
  84.83 93.75 -37.4274.83
Random effects:
  Groups   NameVariance Std.Dev.
  site (Intercept) 0.036592 0.19129
  Residual 1.426886 1.19452
number of obs: 44, groups: site, 15

Fixed effects:
 Estimate Std. Error t value
(Intercept)  -0.4030 0.6857 -0.5877
logpatch  0.1091 0.0491  2.2228
loghab0.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

__
R-help@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.

Re: [R] lmer and estimation of p-values: error with mcmcpvalue()

2007-02-12 Thread Henric Nilsson (Public) 
Den Må, 2007-02-12, 13:58 skrev Christoph Scherber:
 Dear all,

 I am currently analyzing count data from a hierarchical design, and I?ve
 tried to follow the suggestions for a correct estimation of p-values as
 discusssed at R-Wiki
 (http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-testss=lme%20and%20aov).

 However, I have the problem that my model only consists of parameters
 with just 1 d.f. (intercepts, slopes), so that the mcmcpvalue function
 defined below obviously produces error messages.

 How can I proceed in estimating the p-values, then?

 I very much acknowledge any suggestions.

 Best regards
 Christoph.

 ##
 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) }

 m1-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),quasipoisson)

 Generalized linear mixed model fit using Laplace
 Formula: number_pollinators ~ logpatch + loghab + landscape_diversity +
  (1 | site)
 Data: primula
   Family: quasipoisson(log link)
 AIC   BIC logLik deviance
   84.83 93.75 -37.4274.83
 Random effects:
   Groups   NameVariance Std.Dev.
   site (Intercept) 0.036592 0.19129
   Residual 1.426886 1.19452
 number of obs: 44, groups: site, 15

 Fixed effects:
  Estimate Std. Error t value
 (Intercept)  -0.4030 0.6857 -0.5877
 logpatch  0.1091 0.0491  2.2228
 loghab0.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])

Try `mcmcpvalue(as.matrix(markov1[,1]))'.


HTH,
Henric




 Error in colMeans(samp) : 'x' must be an array of at least two dimensions

 __
 R-help@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.



__
R-help@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.

Re: [R] lmer and estimation of p-values: error with mcmcpvalue()

2007-02-12 Thread Christoph Scherber 
Dear Henric,

Thanks, now it works; but how reliable are these estimates? Especially 
with p-values close to 0.05 it is of course important that the range of 
the estimates is not too large. I´ve just run several simulations, each 
of which yielding sometimes quite different p-values.

Best wishes
Christoph


##
Dr. rer. nat. Christoph Scherber
DNPW, Agroecology
University of Goettingen
Waldweg 26
D-37073 Goettingen
http://wwwuser.gwdg.de/~uaoe/Agroecology.html
+49-(0)551-39-8807



Henric Nilsson (Public) schrieb:
 Den Må, 2007-02-12, 13:58 skrev Christoph Scherber:
 Dear all,

 I am currently analyzing count data from a hierarchical design, and I?ve
 tried to follow the suggestions for a correct estimation of p-values as
 discusssed at R-Wiki
 (http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-testss=lme%20and%20aov).

 However, I have the problem that my model only consists of parameters
 with just 1 d.f. (intercepts, slopes), so that the mcmcpvalue function
 defined below obviously produces error messages.

 How can I proceed in estimating the p-values, then?

 I very much acknowledge any suggestions.

 Best regards
 Christoph.

 ##
 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) }

 m1-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),quasipoisson)

 Generalized linear mixed model fit using Laplace
 Formula: number_pollinators ~ logpatch + loghab + landscape_diversity +
  (1 | site)
 Data: primula
   Family: quasipoisson(log link)
 AIC   BIC logLik deviance
   84.83 93.75 -37.4274.83
 Random effects:
   Groups   NameVariance Std.Dev.
   site (Intercept) 0.036592 0.19129
   Residual 1.426886 1.19452
 number of obs: 44, groups: site, 15

 Fixed effects:
  Estimate Std. Error t value
 (Intercept)  -0.4030 0.6857 -0.5877
 logpatch  0.1091 0.0491  2.2228
 loghab0.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])
 
 Try `mcmcpvalue(as.matrix(markov1[,1]))'.
 
 
 HTH,
 Henric
 
 
 
 Error in colMeans(samp) : 'x' must be an array of at least two dimensions

 __
 R-help@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.


 
 
 .


__
R-help@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.

Re: [R] lmer and estimation of p-values: error with mcmcpvalue()

2007-02-12 Thread Gregor Gorjanc 
Christoph Scherber Christoph.Scherber at agr.uni-goettingen.de writes:
 
 Dear Henric,
 
 Thanks, now it works; but how reliable are these estimates? Especially 
 with p-values close to 0.05 it is of course important that the range of 
 the estimates is not too large. I´ve just run several simulations, each 
 of which yielding sometimes quite different p-values.

I have not been using mcmcsamp(), but try increasing number of simulations.
Additionally, since you have only single parameters, you can directly look
at their posterior distributions and compute the areas under the curve.

Gregor

__
R-help@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.