Dear friends, I would like to compare these two non-nested  models:

Alt1 <- glmer(contplt2 ~ 1+A+B+(1|D),family=binomial, data=partNM3)
Alt2 <- glmer(contplt2 ~ 1+A+C+(1|D),family=binomial, data=partNM3)

I tried Faraway(2005) but he only deals with nested models (unless I missed 
something). Can anyone drop a very quick line of help (or just a link which I 
can follow?)

Besides glmer, I tried glmmPQL and MCMCglmm and they all give the same results. 
I would like to thank all the authors of the packages for giving us free 
alternatives. Does anyone know when each of the three is mostly appropriate(or 
they do exactly the same job?)

Thanks

 Dr. Iasonas Lamprianou
Department of Social and Political Sciences
University of Cyprus




________________________________
From: "r-sig-mixed-models-request@r-project.org" 
<r-sig-mixed-models-request@r-project.org>
To: r-sig-mixed-models@r-project.org
Sent: Tue, 24 May, 2011 0:22:09
Subject: R-sig-mixed-models Digest, Vol 53, Issue 25

Send R-sig-mixed-models mailing list submissions to
    r-sig-mixed-models@r-project.org

To subscribe or unsubscribe via the World Wide Web, visit
    https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
    r-sig-mixed-models-request@r-project.org

You can reach the person managing the list at
    r-sig-mixed-models-owner@r-project.org

When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-mixed-models digest..."


Today's Topics:

   1. Re: various priors for MCMCglmm (Ned Dochtermann)
   2. lme/anova vs paired t test query (Ben Bolker)
   3. Re: lme/anova vs paired t test query (Gang Chen)
   4. Re: lme/anova vs paired t test query (Ben Bolker)
   5. Re: lme/anova vs paired t test query (Robert Kushler)


----------------------------------------------------------------------

Message: 1
Date: Mon, 23 May 2011 11:11:07 -0700
From: "Ned Dochtermann" <ned.dochtermann@gmail.com>
To: <r-sig-mixed-models@r-project.org>
Subject: Re: [R-sig-ME] various priors for MCMCglmm
Message-ID: <4ddaa33d.0630440a.2c33.10d2@mx.google.com>
Content-Type: text/plain;    charset="us-ascii"

So in going back over some of this it is clear that my understanding of nu
is still very much lacking, making the multivariate generalizations perhaps
inappropriate. Part of this seems to be that nu incorporates both "degree of
belief" and some aspects of the shape of the prior distribution, at least
for the inverse gamma.

Can anyone recommend some general texts that might clarify some of this?
Gelman and Hill (2006) seem to leave a lot of this out--at least based on my
skimming through it looking for discussion of priors.

Thanks,
Ned




-----Original Message-----
From: Ned Dochtermann [mailto:ned.dochtermann@gmail.com] 
Sent: Friday, May 20, 2011 5:44 PM
To: 'r-sig-mixed-models@r-project.org'
Subject: various priors for MCMCglmm

Hi all,

I have had and continue to have a lot of difficulty with not only the
concept of priors but also their implementation for the MCMCglmm package.
While Jarrod has always been extraordinarily helpful with this--as well as
patient with those of us that seem eternally dense--I thought I'd list out
some of the specifications I've come across to perhaps ease the frequency
with which Jarrod has to reply to questions about priors.

I've specified, for ease, R=G and a single random effect. From things I've
read R is less sensitive to prior specification so this might be okay (?).
Since I'm still not comfortable with priors hopefully these can be clarified
by those more knowledgeable. I'm sure it would also be of general use if
other priors were added via replies.

The parts that are correct should be attributed to Jarrod. The parts that
are wrong should be attributed to my clumsy attempts to generalize something
he wrote. Some of this is just copied from posts to the listserv or Jarrod's
course notes but I thought it might be useful to summarize then in a single
place.

Hopefully none of this is too misleading but will be of some use to others.

Ned

--
Ned Dochtermann
Department of Biology
University of Nevada, Reno

ned.dochtermann@gmail.com
http://wolfweb.unr.edu/homepage/mpeacock/Ned.Dochtermann/
http://www.researcherid.com/rid/A-7146-2010
--



#Inverse Wishart (univariate)

    prior.iw<-list(R=list(V=1, nu=1),G=list(G1=list(V=1, nu=1))) 

#Inverse Wishart (multivariate, variables=x)

    prior.miw<-list(R=list(V=diag(x), nu=x),G=list(G1=list(V=diag(x),
nu=x))) 
#it came up in discussion that this prior can substantially affect the
posterior when there are few data

#I think that the following is an inverse Wishart for multiple response
variables that is less informative is:

    prior.miw<-list(R=list(V=diag(x), nu=x, alpha.mu=c(0,0),
alpha.V=diag(x)*a),G=list(G1=list(V=diag(x), nu=x, alpha.mu=c(0,0),
alpha.V=diag(x)*a)))

#for this a should be "something large (e.g. 1000, depending on the scale of
the data)"

#Inverse Gamma (univariate)

    prior.ig<-list(R=list(V=1, nu=0.002),G=list(G1=list(V=1, nu=0.002)))

#Inverse Gamma (multivariate, two variables)

    prior.mig2<-list(R=list(V=diag(2),
nu=1.002),G=list(G1=list(V=diag(2), nu=1.002))) 

#Inverse Gamma (multivariate, variables=x)

    prior.migX<-list(R=list(V=diag(x),
nu=1.002),G=list(G1=list(V=diag(x), nu=1.002)))

#or 

    prior.migX<-list(R=list(V=diag(x),
nu=b.002),G=list(G1=list(V=diag(x), nu=b.002)))
#where b is x-1

#I initially thought the first of these was the way to go but I'm less so
now after looking at the 
#univariate specification. Also, Gelman (2006. Bayesian Analysis.
1(3):515-533) suggests that
#there can be problems with this prior.

#flat covariance between two responses 

    prior.flatcor<-list(R=list(V=diag(2)*1e-6,
nu=3),G=list(G1=list(V=diag(2)*1e-6, nu=3)))
#to generalize this to x variables V would be diag(x) but I'm not sure if nu
would be x+1

#flat, uniform, prior for just a variance

    prior.flatvar<-list(R=list(V=1e-16, nu=-2),G=list(G1=list(V=1e-16,
nu=-2)))

#another flat improper prior, equivalent to REML fitting:

    prior.reml<-list(G=list(G1=list(V=1,nu=0)),R = list(V =1, nu = ))

#I've actually run the same univariate analysis with prior.reml and lmer and
gotten pretty much the same results



------------------------------

Message: 2
Date: Mon, 23 May 2011 14:55:32 -0400
From: Ben Bolker <bbolker@gmail.com>
To: r-sig-mixed-models <r-sig-mixed-models@r-project.org>
Subject: [R-sig-ME] lme/anova vs paired t test query
Message-ID: <4DDAADA4.4030805@gmail.com>
Content-Type: text/plain; charset=ISO-8859-1

  It was my understanding that testing a two-level fixed effect in a
one-way randomized block ANOVA using REML should be identical to a
paired t test ... but, at least with the attempt below, using
lme(...,method="REML") and t.test(...,paired=TRUE,var.equal=TRUE) seems
to give slightly different results.  They're not different in any
important way, but the difference bothers me.  Can anyone see what I'm
missing?

  Ben Bolker

dat <- structure(list(sex = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("f",
"m"), class = "factor"), prevalence = c(0, 0.375, 0.133333333333333,
0.176470588235294, 0.1875, 0, 0, 1, 1, 0.5, 0.6, 0.333333333333333,
0.5, 0, 0.333333333333333, 0, 0.5, 0, 0.625, 0.333333333333333,
0.5, 0, 0.333333333333333, 0.153846153846154, 0.222222222222222,
0.5, 1, 0.5, 0, 0.277777777777778, 0.125, 0, 0, 0.428571428571429,
0.451612903225806, 0.362068965517241), tripsite = structure(c(1L,
1L, 4L, 4L, 14L, 14L, 5L, 5L, 8L, 8L, 15L, 15L, 6L, 6L, 9L, 9L,
11L, 11L, 16L, 16L, 2L, 2L, 7L, 7L, 10L, 10L, 13L, 13L, 17L,
17L, 3L, 3L, 12L, 12L, 18L, 18L), .Label = c("1.2", "4.2", "5.2",
"1.3", "2.3", "3.3", "4.3", "2.4", "3.4", "4.4", "3.5", "5.5",
"4.6", "1.9", "2.9", "3.9", "4.9", "5.9"), class = "factor")), .Names =
c("sex",
"prevalence", "tripsite"), row.names = c(1L, 2L, 3L, 4L, 9L,
10L, 11L, 12L, 13L, 14L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 38L, 39L, 40L,
41L, 42L, 43L, 45L, 46L), class = "data.frame")


t0 <-
with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))


library(nlme)

(a1 <- anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))

(fstat0 <- t0$statistic^2)     ## 0.789627
(fstat1 <- a1[["F-value"]][2]) ## 0.8056624

tmpf <- function(dat) {
  t0 <-
with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))

  (a1 <-
anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
  (fstat0 <- t0$statistic^2)     ## 0.789627
  (fstat1 <- a1[["F-value"]][2]) ## 0.8056624
  fstat0-fstat1
}

  If I then redo this with random data, the results are within 1e-5
set.seed(1001)
rr <- replicate(500,tmpf(transform(dat,prevalence=rnorm(nrow(dat)))))
##
##  summary(rr)
##      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
## -1.7770000 -0.0568400 -0.0000015 -0.0722900  0.0000000  0.0000139

mean(abs(rr)<1e-5)
## 0.502

hist(log10(abs(rr)),breaks=50,col="gray")

## this is sort of interesting, suggesting a bimodal distribution with
one component (the 'good ones' centered at 10^{-8} and another centered
at 0.1 ...



------------------------------

Message: 3
Date: Mon, 23 May 2011 15:40:29 -0400
From: Gang Chen <gangchen@mail.nih.gov>
To: Ben Bolker <bbolker@gmail.com>
Cc: r-sig-mixed-models <r-sig-mixed-models@r-project.org>
Subject: Re: [R-sig-ME] lme/anova vs paired t test query
Message-ID: <BANLkTikDd0bN2Bsq6G31VHsNXV=t0qw11g@mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1

I believe I know the answer for this one because I used to struggle
with this for quite a while. Eventually I figured it out.

To match up with the paired t-test, you need

(a2 <- 
anova(lme(prevalence~sex,random=list(tripsite=pdCompSymm(~sex-1)),data=dat,method="REML")))

(fstat2 <- a2[["F-value"]][2]) # you get 0.789627

The catch here is that the compound-symmetry structure for the
variance-covariance structure of the random components in the lme
model is equivalent to a paired t-test, a special of the traditional
repeated-meaures one-way ANOVA.

Gang


On Mon, May 23, 2011 at 2:55 PM, Ben Bolker <bbolker@gmail.com> wrote:
> ?It was my understanding that testing a two-level fixed effect in a
> one-way randomized block ANOVA using REML should be identical to a
> paired t test ... but, at least with the attempt below, using
> lme(...,method="REML") and t.test(...,paired=TRUE,var.equal=TRUE) seems
> to give slightly different results. ?They're not different in any
> important way, but the difference bothers me. ?Can anyone see what I'm
> missing?
>
> ?Ben Bolker
>
> dat <- structure(list(sex = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("f",
> "m"), class = "factor"), prevalence = c(0, 0.375, 0.133333333333333,
> 0.176470588235294, 0.1875, 0, 0, 1, 1, 0.5, 0.6, 0.333333333333333,
> 0.5, 0, 0.333333333333333, 0, 0.5, 0, 0.625, 0.333333333333333,
> 0.5, 0, 0.333333333333333, 0.153846153846154, 0.222222222222222,
> 0.5, 1, 0.5, 0, 0.277777777777778, 0.125, 0, 0, 0.428571428571429,
> 0.451612903225806, 0.362068965517241), tripsite = structure(c(1L,
> 1L, 4L, 4L, 14L, 14L, 5L, 5L, 8L, 8L, 15L, 15L, 6L, 6L, 9L, 9L,
> 11L, 11L, 16L, 16L, 2L, 2L, 7L, 7L, 10L, 10L, 13L, 13L, 17L,
> 17L, 3L, 3L, 12L, 12L, 18L, 18L), .Label = c("1.2", "4.2", "5.2",
> "1.3", "2.3", "3.3", "4.3", "2.4", "3.4", "4.4", "3.5", "5.5",
> "4.6", "1.9", "2.9", "3.9", "4.9", "5.9"), class = "factor")), .Names =
> c("sex",
> "prevalence", "tripsite"), row.names = c(1L, 2L, 3L, 4L, 9L,
> 10L, 11L, 12L, 13L, 14L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
> 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 38L, 39L, 40L,
> 41L, 42L, 43L, 45L, 46L), class = "data.frame")
>
>
> t0 <-
>with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
>)
>
> library(nlme)
>
> (a1 <- anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
>
> (fstat0 <- t0$statistic^2) ? ? ## 0.789627
> (fstat1 <- a1[["F-value"]][2]) ## 0.8056624
>
> tmpf <- function(dat) {
> ?t0 <-
>with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
>)
> ?(a1 <-
> anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
> ?(fstat0 <- t0$statistic^2) ? ? ## 0.789627
> ?(fstat1 <- a1[["F-value"]][2]) ## 0.8056624
> ?fstat0-fstat1
> }
>
> ?If I then redo this with random data, the results are within 1e-5
> set.seed(1001)
> rr <- replicate(500,tmpf(transform(dat,prevalence=rnorm(nrow(dat)))))
> ##
> ## ?summary(rr)
> ## ? ? ?Min. ? ?1st Qu. ? ? Median ? ? ? Mean ? ?3rd Qu. ? ? ? Max.
> ## -1.7770000 -0.0568400 -0.0000015 -0.0722900 ?0.0000000 ?0.0000139
>
> mean(abs(rr)<1e-5)
> ## 0.502
>
> ?hist(log10(abs(rr)),breaks=50,col="gray")
>
> ## this is sort of interesting, suggesting a bimodal distribution with
> one component (the 'good ones' centered at 10^{-8} and another centered
> at 0.1 ...
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



------------------------------

Message: 4
Date: Mon, 23 May 2011 16:10:16 -0400
From: Ben Bolker <bbolker@gmail.com>
To: Gang Chen <gangchen@mail.nih.gov>
Cc: r-sig-mixed-models <r-sig-mixed-models@r-project.org>
Subject: Re: [R-sig-ME] lme/anova vs paired t test query
Message-ID: <4DDABF28.8020504@gmail.com>
Content-Type: text/plain; charset=ISO-8859-1

  Great, very clear!  Thanks.

  Ben

On 05/23/2011 03:40 PM, Gang Chen wrote:
> I believe I know the answer for this one because I used to struggle
> with this for quite a while. Eventually I figured it out.
> 
> To match up with the paired t-test, you need
> 
> (a2 <- 
>anova(lme(prevalence~sex,random=list(tripsite=pdCompSymm(~sex-1)),data=dat,method="REML")))
>
> (fstat2 <- a2[["F-value"]][2]) # you get 0.789627
> 
> The catch here is that the compound-symmetry structure for the
> variance-covariance structure of the random components in the lme
> model is equivalent to a paired t-test, a special of the traditional
> repeated-meaures one-way ANOVA.
> 
> Gang
> 
> 
> On Mon, May 23, 2011 at 2:55 PM, Ben Bolker <bbolker@gmail.com> wrote:
>>  It was my understanding that testing a two-level fixed effect in a
>> one-way randomized block ANOVA using REML should be identical to a
>> paired t test ... but, at least with the attempt below, using
>> lme(...,method="REML") and t.test(...,paired=TRUE,var.equal=TRUE) seems
>> to give slightly different results.  They're not different in any
>> important way, but the difference bothers me.  Can anyone see what I'm
>> missing?
>>
>>  Ben Bolker
>>
>> dat <- structure(list(sex = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
>> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
>> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("f",
>> "m"), class = "factor"), prevalence = c(0, 0.375, 0.133333333333333,
>> 0.176470588235294, 0.1875, 0, 0, 1, 1, 0.5, 0.6, 0.333333333333333,
>> 0.5, 0, 0.333333333333333, 0, 0.5, 0, 0.625, 0.333333333333333,
>> 0.5, 0, 0.333333333333333, 0.153846153846154, 0.222222222222222,
>> 0.5, 1, 0.5, 0, 0.277777777777778, 0.125, 0, 0, 0.428571428571429,
>> 0.451612903225806, 0.362068965517241), tripsite = structure(c(1L,
>> 1L, 4L, 4L, 14L, 14L, 5L, 5L, 8L, 8L, 15L, 15L, 6L, 6L, 9L, 9L,
>> 11L, 11L, 16L, 16L, 2L, 2L, 7L, 7L, 10L, 10L, 13L, 13L, 17L,
>> 17L, 3L, 3L, 12L, 12L, 18L, 18L), .Label = c("1.2", "4.2", "5.2",
>> "1.3", "2.3", "3.3", "4.3", "2.4", "3.4", "4.4", "3.5", "5.5",
>> "4.6", "1.9", "2.9", "3.9", "4.9", "5.9"), class = "factor")), .Names =
>> c("sex",
>> "prevalence", "tripsite"), row.names = c(1L, 2L, 3L, 4L, 9L,
>> 10L, 11L, 12L, 13L, 14L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
>> 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 38L, 39L, 40L,
>> 41L, 42L, 43L, 45L, 46L), class = "data.frame")
>>
>>
>> t0 <-
>>with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
>>)
>>
>> library(nlme)
>>
>> (a1 <- anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
>>
>> (fstat0 <- t0$statistic^2)     ## 0.789627
>> (fstat1 <- a1[["F-value"]][2]) ## 0.8056624
>>
>> tmpf <- function(dat) {
>>  t0 <-
>>with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
>>)
>>  (a1 <-
>> anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
>>  (fstat0 <- t0$statistic^2)     ## 0.789627
>>  (fstat1 <- a1[["F-value"]][2]) ## 0.8056624
>>  fstat0-fstat1
>> }
>>
>>  If I then redo this with random data, the results are within 1e-5
>> set.seed(1001)
>> rr <- replicate(500,tmpf(transform(dat,prevalence=rnorm(nrow(dat)))))
>> ##
>> ##  summary(rr)
>> ##      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
>> ## -1.7770000 -0.0568400 -0.0000015 -0.0722900  0.0000000  0.0000139
>>
>> mean(abs(rr)<1e-5)
>> ## 0.502
>>
>>  hist(log10(abs(rr)),breaks=50,col="gray")
>>
>> ## this is sort of interesting, suggesting a bimodal distribution with
>> one component (the 'good ones' centered at 10^{-8} and another centered
>> at 0.1 ...
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>



------------------------------

Message: 5
Date: Mon, 23 May 2011 17:21:54 -0400
From: Robert Kushler <kushler@oakland.edu>
To: r-sig-mixed-models@r-project.org
Subject: Re: [R-sig-ME] lme/anova vs paired t test query
Message-ID: <4DDACFF2.7090103@oakland.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed

Ben,

The "good ones" in your simulation correspond to cases where the classical
"ANOVA" estimate of the block variance component is greater than zero (i.e.
the F statistic for "tripsite" in the two way fixed effects anova is greater
than 1).

> tmpf <- function(dat) {
+   t0 <- 
with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))

+   a1 <- anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML"))
+   a2 <- anova(lm(prevalence~sex+tripsite,data=dat))
+   fstat0 <- t0$statistic^2     ## 0.789627
+   fstat1 <- a1[["F-value"]][2] ## 0.8056624
+   fstat2 <- a2[["F value"]][2]
+   c(fstat0,fstat1,fstat2)
+ }
>
> set.seed(1001)
> rr <- unlist(replicate(500,tmpf(transform(dat,prevalence=rnorm(nrow(dat))))))
> table((abs(rr[1,]-rr[2,])<0.00002),(rr[3,]>=1))

         FALSE TRUE
   FALSE   246    0
   TRUE      3  251



Regards,   Rob Kushler


On 5/23/2011 4:10 PM, Ben Bolker wrote:
>    Great, very clear!  Thanks.
>
>    Ben
>
> On 05/23/2011 03:40 PM, Gang Chen wrote:
>> I believe I know the answer for this one because I used to struggle
>> with this for quite a while. Eventually I figured it out.
>>
>> To match up with the paired t-test, you need
>>
>> (a2<- 
>>anova(lme(prevalence~sex,random=list(tripsite=pdCompSymm(~sex-1)),data=dat,method="REML")))
>>
>> (fstat2<- a2[["F-value"]][2]) # you get 0.789627
>>
>> The catch here is that the compound-symmetry structure for the
>> variance-covariance structure of the random components in the lme
>> model is equivalent to a paired t-test, a special of the traditional
>> repeated-meaures one-way ANOVA.
>>
>> Gang
>>
>>
>> On Mon, May 23, 2011 at 2:55 PM, Ben Bolker<bbolker@gmail.com>  wrote:
>>>   It was my understanding that testing a two-level fixed effect in a
>>> one-way randomized block ANOVA using REML should be identical to a
>>> paired t test ... but, at least with the attempt below, using
>>> lme(...,method="REML") and t.test(...,paired=TRUE,var.equal=TRUE) seems
>>> to give slightly different results.  They're not different in any
>>> important way, but the difference bothers me.  Can anyone see what I'm
>>> missing?
>>>
>>>   Ben Bolker
>>>
>>> dat<- structure(list(sex = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
>>> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
>>> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("f",
>>> "m"), class = "factor"), prevalence = c(0, 0.375, 0.133333333333333,
>>> 0.176470588235294, 0.1875, 0, 0, 1, 1, 0.5, 0.6, 0.333333333333333,
>>> 0.5, 0, 0.333333333333333, 0, 0.5, 0, 0.625, 0.333333333333333,
>>> 0.5, 0, 0.333333333333333, 0.153846153846154, 0.222222222222222,
>>> 0.5, 1, 0.5, 0, 0.277777777777778, 0.125, 0, 0, 0.428571428571429,
>>> 0.451612903225806, 0.362068965517241), tripsite = structure(c(1L,
>>> 1L, 4L, 4L, 14L, 14L, 5L, 5L, 8L, 8L, 15L, 15L, 6L, 6L, 9L, 9L,
>>> 11L, 11L, 16L, 16L, 2L, 2L, 7L, 7L, 10L, 10L, 13L, 13L, 17L,
>>> 17L, 3L, 3L, 12L, 12L, 18L, 18L), .Label = c("1.2", "4.2", "5.2",
>>> "1.3", "2.3", "3.3", "4.3", "2.4", "3.4", "4.4", "3.5", "5.5",
>>> "4.6", "1.9", "2.9", "3.9", "4.9", "5.9"), class = "factor")), .Names =
>>> c("sex",
>>> "prevalence", "tripsite"), row.names = c(1L, 2L, 3L, 4L, 9L,
>>> 10L, 11L, 12L, 13L, 14L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
>>> 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 38L, 39L, 40L,
>>> 41L, 42L, 43L, 45L, 46L), class = "data.frame")
>>>
>>>
>>> t0<-
>>>with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
>>>)
>>>
>>> library(nlme)
>>>
>>> (a1<- anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
>>>
>>> (fstat0<- t0$statistic^2)     ## 0.789627
>>> (fstat1<- a1[["F-value"]][2]) ## 0.8056624
>>>
>>> tmpf<- function(dat) {
>>>   t0<-
>>>with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
>>>)
>>>   (a1<-
>>> anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML")))
>>>   (fstat0<- t0$statistic^2)     ## 0.789627
>>>   (fstat1<- a1[["F-value"]][2]) ## 0.8056624
>>>   fstat0-fstat1
>>> }
>>>
>>>   If I then redo this with random data, the results are within 1e-5
>>> set.seed(1001)
>>> rr<- replicate(500,tmpf(transform(dat,prevalence=rnorm(nrow(dat)))))
>>> ##
>>> ##  summary(rr)
>>> ##      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
>>> ## -1.7770000 -0.0568400 -0.0000015 -0.0722900  0.0000000  0.0000139
>>>
>>> mean(abs(rr)<1e-5)
>>> ## 0.502
>>>
>>>   hist(log10(abs(rr)),breaks=50,col="gray")
>>>
>>> ## this is sort of interesting, suggesting a bimodal distribution with
>>> one component (the 'good ones' centered at 10^{-8} and another centered
>>> at 0.1 ...
>>>
>>> _______________________________________________
>>> R-sig-mixed-models@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



------------------------------

_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


End of R-sig-mixed-models Digest, Vol 53, Issue 25
**************************************************

	[[alternative HTML version deleted]]


