# Identifying subject specific outliers in the presence of between subject heterogeneity

I have a simple dataset of people's heights, many people with measurements on multiple days (once a year for 10 years, say). I have the date of each measurement.

Some of the height values are absurd. I already drop values that are 'impossible' (e.g. height values above 3m). However, I'd like to identify unusual values within a patient. If a patient has 5 records at around 1.8m, I'd like a red flag if there's a value for 2.1m

I am considering trying the ESD test for outliers - it seems straightforward to implement - but I thought I'd ask if there are better ideas before I get started.

Thanks.

EDIT:

The ESD test I refer to is the [[Generalized] 'Extreme Studentized Deviate'][1] test; similar to the Grubbs Test, it 'is used to detect one or more outliers in a univariate data set that follows an approximately normal distribution'.

I'll have measurements for over 100,000 people - however, each person will only have between 0 and around 20 measurements.

I'm hoping to generalise this to other measurements such as BMI, lab tests (Eosinophils, White Blood Counts), respiratory tests (FEV1), etc. I suspect the latter ones will be difficult as they're more varied, rising and falling with little dependency on time - hence my starting with height!

• What is the "ESD test"? Aug 29, 2014 at 15:50
• what is the sample size --the number of people you have measurements for? Aug 29, 2014 at 18:45
• I've edited the original comment to describe the ESD test. Sep 1, 2014 at 7:29

I think you will have to take a look at random effects (see How to account for participants in a study design?).

The idea is that you have, for a number of people $N$, a number of height measurements in different years $y$, so $h_{iy},i=1, 2 \dots, N, y=1, 2, \dots 10$. Each $h_{iy}$ depends on the height in the first year (which depends on the person $i$) and increases with the years (the rate of increase can also depend on the person $i$), so you have something like $h_{iy} = \beta_{1i} + \beta_{2i} y + \epsilon$, where the coefficients $\beta_{1i}, \beta_{2i}$ depend on the person (hence the subscipt $i$), but if you can estimate such an equation, then you have one regression line for each person and you can look, for that single person, which points are 'far away' from that individual person's own regression line.

A reference book is Fitzmaurice, Laird, Ware, "Applied longitudinal analysis" where you can find all the details, so I will just explain the main lines of the reasoning:

It is assumed that the $\beta_{1i}$ and $\beta_{2i}$ are normally distributed with a mean $\beta_1$ (resp. $\beta_2$) and a standard deviation $\sigma_1$ (resp. $\sigma_2$) (so the assumption is that the heights in the first year are normally distributed, just as the growth rates). So we can re-write $\beta_{ki}$ as $\beta_{ki}=\beta_k + b_{kj}$ where $\beta_k$ is the mean of all $\beta_{ki}$ and $b_{ki} \sim N(0;\sigma_i), k=1,2$.

So your regression becomes $h_{iy}=\beta_1+b_{1i} + (\beta_2+b_{2i})y+\epsilon$, this is the form that you find in the reference Fitzmaurice et al.

the $b_{ki}$ are called the random effects. The $\beta_{k}$ are the fixed effects.

To illustrate this let me use some simulated data:

# this is just for simuating data, you have the data, the simulation generates heights
# for N=500 people for 10 years.

library(reshape2)

# simulation for 500 persons
N<-500

#simulate the heights in the first year, mean 1.7, sd 0.2,
# generate a random number of each of the N persons
h0<-rnorm(n=N, mean=1.7, sd=0.2)

#generate a random growth rate for each of the N persons
ch1<-rnorm(n=N, mean=0.01, sd=0.01)

# simulate h_iy:  the height in the first year + (annual growth x year) + error term
h<-h0 + ch1 %*% t(1:10) + rnorm(n=N, mean=0, sd=0.02)

# reformat it a bit so that it is useable later on
colnames(h)<-1:10
df<-cbind( data.frame(subject=1:nrow(h)),
as.data.frame(h) )

df.molten<-melt(data=df, id.vars="subject", value.name="height", variable.name="year")
df.molten$year<-as.numeric(df.molten$year)


After simulating the data, we will now estimate the above mentioned equation as a linear mixed effects model (with the nlme package in R):

library(nlme)
library(ggplot2)

lme.h<-lme(height  ~ year + 1,       # we estimate height as function of year
data=df.molten,           # on the simulated data
random= ~ 1 + year | subject,  # we use a random effect on the intercept and on the coeffcient of year
control=lmeControl(opt="optim"),
method="REML")


The values that are estimated for your data can be found with the following calls:

lme.fix<-fixed.effects(lme.h)
lme.blup<-random.effects(lme.h)


The first one gives you the fixed effects $\beta_k$ (compare them to the values used in the simulation), the second one gives you the best linear unbiased prediction (blup) of the random effect for each person of your dataset.

Compare now the 'estimated fixed effects + the blup predictors of the random effects' for each person to your input (in this case simulated) data, first for the intercept:

df.plot<-data.frame(sim=h0, estim=lme.blup[[1]]+lme.fix[1])

ggplot(df.plot, aes(x=sim,y=estim))+
coord_fixed()+
geom_abline(intercept=0, slope=1, colour="red")+
geom_point()


I can not upload plots from here but if you execute the code then you get a graph with on the horizontal axis the simulated initial (i.e. in the first year) heights and on the vertical axis the estimated/predicted values for each person in the first year, the red line is the bisector.

Similar, but for the coefficient of the year:

df.plot<-data.frame(sim=ch1, estim=lme.blup[[2]]+lme.fix[2])

ggplot(df.plot, aes(x=sim,y=estim))+
coord_fixed()+
geom_abline(intercept=0, slope=1, colour="red")+
geom_point()


With these results you have a regression line for each inidivual person, you can also estimate the variance around the line and then find ''outliers'' for each person's line. For details I refer to Fitzmaurice et al.

so for person $i$ you will find the regression line $h_{iy}=\hat{\beta}_1+\hat{b}_{1i} + (\hat{\beta}_2+\hat{b}_{2i})y$ where $\hat{\beta}_1$ is in lme.fix[1], $\hat{b}_{1i}$ in lme.blup[[1]][i] and $\hat{\beta}_2$ is in lme.fix[2], $\hat{b}_{2i}$ in lme.blup[[2]][i]

For the details on the nlme package you can use this pdf

Since you have a measurement that is very well known, why not define deviation based on your knowledge of measurement technology rather than use a statistical measure of deviance. For example, if you have experience making height measurements and you know height easily can be measured accurately to within 1" (about .05 m), then you might simply discard all measures more deviant than 1" as flawed. Statistical analysis is about finding meaning in data. If you have a definition of a deviant measure based on logical evaluation of the limits of height measurement technology, why look to statistical analysis for a measure of deviance?

• the measures are taken over 10 years Sep 1, 2014 at 18:06
• Likely there is an existing, published literature on weight change over time. Your definition of a large (suspect) change could be based on that literature. Sep 1, 2014 at 18:36
• I would be very cautious using published literature on variable change over time to define large change. Your subpopulation will, necessarily, be different than the the other study's subpopulation and it may be reasonable that this difference affects rates of change: for example, if you are looking at a novel diet, you want not want to identify unusual change from estimates coming from a different study. Not saying that there are no measurements for which this approach could work, but it should be used with some caution. Aug 16, 2015 at 20:39