Skip to content

Question concerning anova()

6 messages · Bin.You at Forst.bwl.de, Bert Gunter, Rainer M Krug

Rainer M Krug
# 
Hi

I am comparing four different linear mixed effect models, derived from updating the original one. To 
compare these, I want to use anova(). I therefore do the following (not reproducible - just to 
illustration purpose!):

dat <- loadSPECIES(SPECIES)
subs <- expression(dead==FALSE & recTreat==FALSE)
feff <- noBefore~pHarv*year      # fixed effect in the model
reff <- ~year|plant              # random effect in the model, where year is the
corr <- corAR1(form=~year|plant) # describing the within-group correlation structure
#
dat.lme <- lme(
              fixed = feff,                           # fixed effect in the model
              data  = dat,
              subset = eval(subs),
              method = "ML",
              random = reff,                          # random effect in the model
              correlation = corr,
              na.action = na.omit
              )
dat.lme.r1 <- update(dat.lme, random=~1|plant)
dat.lme.f1 <- update(dat.lme, fixed=noBefore~year)
dat.lme.r1.f1 <- update(dat.lme.r1, fixed=noBefore~year)


The anova is as follow:

 > anova(dat.lme, dat.lme.r1, dat.lme.f1, dat.lme.r1.f1)
               Model df      AIC      BIC    logLik   Test      L.Ratio p-value
dat.lme           1  9 1703.218 1733.719 -842.6089
dat.lme.r1        2  7 1699.218 1722.941 -842.6089 1 vs 2 1.019230e-07       1
dat.lme.f1        3  7 1705.556 1729.279 -845.7779
dat.lme.r1.f1     4  5 1701.556 1718.501 -845.7779 3 vs 4 8.498318e-08       1

I have two questions:
1) I am wondering why the "2 vs 3" does not give the Test values?
Is this because the two models are considered as "identical", which would be strange, due to the 
different logLik values.

2) If I want to compare all models among each other - is there a "best" way? I would be reluctant to 
do several ANOVA's, due to necessary corrections for multple tests (although this should not be a 
problem here?)

I can obviously select the best model based on the AIC.

Thanks in advance,

Rainer
Bert Gunter
# 
Models with different fixed effects estimated by REML cannot be
compared by anova.

In future, please post questions on mixed effects models on the
r-sig-mixed-effects mailing lists. You're likely to receive more
informative replies there, too.

-- Bert
On Wed, Aug 22, 2012 at 7:23 AM, Rainer M Krug <r.m.krug at gmail.com> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Bin.You at Forst.bwl.de
    

    

      
      #
    

  


  

      
Hi, 

I am doing Bayesian Calibration for a model which has 24 parameters using MCMC. I get a sample of 100,000 points from the posterior pdf. Let's say, It is a dataframe named as "BC_output", with 24 parameters (V1,V2,....V24)

I would like to show the Maximum a Posteriori (MAP)parameter value of the model. 

Can I do it like this: 

For V1: density(BC_output$V1)$x[which.max(density(BC_output$V1)$y)]

If not, what is the correct function to estimate a MAP value in R? 

Thanks a lot,

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rainer M Krug
    

    

      
      #
    

  


  

      
On 22/08/12 16:36, Bert Gunter wrote:

            

      
      
      
    

        
I have seen that much in "Modern Applied Statistics in S", and therefore have chosen the model = "ML"

            

      
      
      
    

        
Thanks - wasn't aware of this sig - I'll send the reply there as well.

Thanks,

Rainer

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Bert Gunter
    

    

      
      #
    

  


  

      
Oops -- missed that. OTOH, my reply demonstrates the value of the
mixed models list recommendation.

-- Bert

        
On Wed, Aug 22, 2012 at 7:55 AM, Rainer M Krug <r.m.krug at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rainer M Krug
    

    

      
      #
    

  


  

      
Further discussed on r-sig-mixed-models

Rainer

        
On 22/08/12 17:04, Bert Gunter wrote: