Skip to content

Calculating repeatability with ordinal data

7 messages · Federico Calboli, Jarrod Hadfield, Pierre de Villemereuil +1 more

Original
Samantha Patrick
# 
Hi

Firstly I sent this message last week but can find no evidence that it actually sent.
However, if this is a double posting I am very sorry.

I am currently working with repeated measures for individuals and I am
trying to quantify individual repeatability.  Normally, for continuous
distributions, I would use a mixed model and calculate the variance
explained by individual divided by the total variance.  However my individual scores
are ordinal and I have been using the clmm function in the Ordinal package:

example of data:
ID        Bird        Sex    scoremax    year
622    BS8831    M       2                 2008
623    BS8831    M       1                 2010
624    BS8831    M       1                 2011
625    BS9065    M       1                 2010
626    BS9065    M       3                 2011
627   BS19724    F       4                 2010
628   BS19724    F       5                 2010
629   BS21302    F       1                 2010
630   BS25376    F       1                 2011
631    BS9184    F       2                 2009
632   BS19989    M       3                 2011
633   BS21617    M       4                 2008
634   BS21617    M       2                 2009
635   BS21617    M       1                 2010

where scoremax ranges from 1-5, and there are 1188 birds and 1638
observations.

scoremax<-as.factor(scoremax)
bird<-as.factor(bird)
year<-as.factor(year)
fmm1<- clmm(scoremax~year+ (1|bird), link = c("probit"), Hess =TRUE)

summary(fmm1)

but this only gives the variance estimate for bird, with no residual
estimate.  Some investigations reveal that using an ordinal regression
in MCMCglmm will also not estimate the residual variance, and it seems
you need to constrain this value.  I have been unable to find any posts
about repeatability in ordinal data.

My questions I guess are:
Is using a mixed model appropriate for calculating the repeatability of
ordinal data (and if not does anyone know any other methods)?

If it is, does anyone have any hints on how to calculate the residual variance,
to enable repeatability estimates to be calculated.

Many Thanks

Sam
Federico Calboli
#
On 2 May 2012, at 11:08, Samantha Patrick wrote:

            

      
      
      
    

        
does 

scoremax = as.ordered(scoremax)

make any difference in your results?

I ask this because as.factor() does not, strictly speaking, create an ordered factor.

BW

F

            

      
      
      
    

        

      
      
    

        
--
Federico C. F. Calboli
Neuroepidemiology and Ageing Research
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

  
  


  


      

    

    

      
      

        


  

        

      Samantha Patrick
    

    

      
      #
    

  


  

      
Hi

Thank you for the tip - I had editted my code!  However, sadly it 
doesn't  solve the problem of how to estimate the residual variance.

Many Thanks

Sam


Le 02/05/2012 12:38, Federico Calboli a ?crit :

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi,

The residual variance can never be estimated from ordinal data. Most  
programs will set it to zero, some programs allow you to set it at  
anything (MCMCglmm). I have not seen repeatabilities for ordinal data  
but I presume you can add the variance of the relevant distribution  
(pi^2/3 or 1 for logit/probit) to the denominator in order to get an  
intra-class correlation. I'm not confident how this would be  
interpreted for an ordinal response though. Note that if you do use  
MCMCglmm you need to include the constrained residual variance in the  
denominator (i.e. to get the denominator you need to add two to the  
"Bird" variance if you have constrained the residual variance to one).

Cheers,

Jarrod


Quoting Samantha Patrick <spatrick at cebc.cnrs.fr> on Wed, 02 May 2012  
13:11:09 +0200:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
20 days later

    

      
      

        


  

        

      Samantha Patrick
    

    

      
      #
    

  


  

      
Hi

Thanks for the comments - sorry for the delayed response, I had to learn 
how to use MCMCglmm!

Following Shinichi's advice I have been using an adapted version of the 
code in Dean et al 2012 (Am Nat). To start with I have been trying to 
fit the data as categorical, as the original script is based on this 
distribution.

Simplified code:

Bird = factor with 1183 levels with 1638 observations
Year = factor with 4 levels
scoremax = ordinal factor with 5 levels (0-4)


year<-as.factor(year)

# setting prior probabilities for random effects
IJ<-(1/5)*(diag(4) + matrix(1,4,4))

priortest1 <- list(R=list(V = IJ, nu = 0.002), G = list(G1 = list(V = 1, 
nu = 0.002), G2 = list(V = 1, nu = 0.002)))

model1b<-MCMCglmm(scoremax~trait-1, random =~Bird +year, rcov = 
~us(trait):units, data = Data, prior = priortest1, family = "categorical")

# Getting "adjusted" repeatability - Equation A10
summary(model1$VCV[,1]/(model1$VCV[,1]+2/3+pi^2/3))
posterior.mode(model1$VCV[,1]/(model1$VCV[,1]+2/3+pi^2/3))

This does run but the model fit is pretty bad.  I haven't tried to 
optimise it as it is important that my data is ordered so I am trying to 
work on the script below:

priortest2 <- list(R=list(V = 1 , fix=1), G = list(G1=list(V=1, nu=1, 
alpha.mu=0, alpha.V=100),G2=list(V=1, nu=1, alpha.mu=0, alpha.V=100))) 
##  this seems to be a commonly used prior for ordinal data so I was 
using it as a starting point but have not got as far as to test whether 
it is more suitable than priortest1 due to the problem below

model1b<-MCMCglmm(scoremax~trait-1, random =~Bird + year, rcov = 
~us(trait):units, data = Data, prior = priortest2, family = "ordinal")

At this stage I get an error message which says:

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
   contrasts can be applied only to factors with 2 or more levels

I have tried to create dummy data to seek out the problem.  I thought 
the problem may stem from the model trying to estimate contrasts between 
each level of the ordinal response ie. 0-1, 1-0, 0-2, 2-0, 0-3, 3-0, 0-4 
etc.  Within the data individuals move between all scores but not both 
ways (i.e.  There are individuals who change between 1 and 2, but not 
necessarily individuals that move from 1 to 2 AND 2 to 1).  However, I 
have not fitted observation number in the model (the temporal order of 
observations is not specified) and altering the order of observation 
does not seem to fix this.  I have simplified the model and run without 
year, but I get the same error.

I can happily post more data if it is helpful.

Any advice would be greatly appreciated!

Thanks

Sam




Le 02/05/2012 15:43, Jarrod Hadfield a ?crit :