title: Identifying outlying datapoints from residuals (GeoLight package)
body: <p>I am analysing some data collected from a geolocator placed on a migratory bird. In a nutshell, my data are sunrise and sunset times, which are then used to determine position on the globe.</p>
  
  <p>I am using a package GeoLight (<a href="http://cran.r-project.org/web/packages/GeoLight/GeoLight.pdf" rel="nofollow">http://cran.r-project.org/web/packages/GeoLight/GeoLight.pdf</a>) to identify outlying data - specifically, I am using the LoessFilter function which applies a polynomial regression and identify residuals that are greater than 3 interquantile ranges (specified by k in the code when applying the function)</p>
  
  <p><strong>My problem is</strong>: the function returns graphs in which outlying datapoints are identified in red. There seems to be an issue with the code itself regarding returned TRUE or FALSE statements stating which points are outliers - all are stated as TRUE, even if outliers are identified.</p>
  
  <p>I have therefore modified the function code to state which residuals are outliers.</p>
  
  <p>However, when I then remove those rows from the original dataset and re-run the function, the points have not been removed. Therefore, there is some discrepancy between which residuals are relating to values in the original data: i.e. if the output states that residual 78 is an outlying point, removing row 78 from the original data does not remove the outlying datapoint.</p>
  
  <p>I would very much appreciate some help with removing the outlying datapoints identified using the function. It seems like a very easy fix but I can't seem to figure it out.</p>
  
  <p>Code for full function and data below</p>
  
  <p>Thanks
  Emma</p>
  
  <pre><code>log2$tFirst&lt;-as.POSIXlt(log2$tFirst)
  log2$tSecond&lt;-as.POSIXlt(log2$tSecond)
  </code></pre>
  
  <p>CODE TO GET OUTLYING RESIDUALS</p>
  
  <pre><code>i.get.outliers&lt;-function(residuals, k=3) {
    x &lt;- residuals
  </code></pre>
  
  <p># x is a vector of residuals
    # k is a measure of how many interquartile ranges to take before saying that point is an outlier
    # it looks like 3 is a good preset for k</p>
  
  <pre><code>QR&lt;-quantile(x, probs = c(0.25, 0.75))
  IQR&lt;-QR[2]-QR[1]
  Lower.band&lt;-QR[1]-(k*IQR)
  Upper.Band&lt;-QR[2]+(k*IQR)
  delete&lt;-which(x&lt;Lower.band |  x&gt;Upper.Band)
  return(as.vector(delete))
  }
  </code></pre>
  
  LOESS FILTER FUNCTION CODE
  
  <pre><code>loessFilter &lt;- function(tFirst, tSecond, type, k=3, plot=TRUE){
  tw &lt;- data.frame(datetime=as.POSIXct(c(tFirst,tSecond),"UTC"),type=c(type,ifelse(type==1,2,1)))
  tw &lt;- tw[!duplicated(tw$datetime),]
  tw &lt;- tw[order(tw[,1]),]
  
  hours &lt;- as.numeric(format(tw[,1],"%H"))+as.numeric(format(tw[,1],"%M"))/60
  
  for(t in 1:2){
  cor &lt;- rep(NA, 24)
  for(i in 0:23){
    cor[i+1] &lt;- max(abs((c(hours[tw$type==t][1],hours[tw$type==t])+i)%%24 -
                          (c(hours[tw$type==t],hours[tw$type==t][length(hours)])+i)%%24),na.rm=T)
  }
  hours[tw$type==t] &lt;- (hours[tw$type==t] + (which.min(round(cor,2)))-1)%%24
  }
  
  dawn &lt;- data.frame(id=1:sum(tw$type==1),
                   datetime=tw$datetime[tw$type==1],
                   type=tw$type[tw$type==1],
                   hours = hours[tw$type==1], filter=FALSE)
  
  dusk &lt;- data.frame(id=1:sum(tw$type==2),
                   datetime=tw$datetime[tw$type==2],
                   type=tw$type[tw$type==2],
                   hours = hours[tw$type==2], filter=FALSE)
  
  
  for(d in seq(30,k,length=5)){
  
  predict.dawn &lt;- predict(loess(dawn$hours[!dawn$filter]~as.numeric(dawn$datetime[!dawn$filter]),span=0.1))
  predict.dusk &lt;- predict(loess(dusk$hours[!dusk$filter]~as.numeric(dusk$datetime[!dusk$filter]),span=0.1))
  
  del.dawn &lt;- i.get.outliers(as.vector(residuals(loess(dawn$hours[!dawn$filter]~
                                                         as.numeric(dawn$datetime[!dawn$filter]),span=0.1))),k=d)
  del.dusk &lt;- i.get.outliers(as.vector(residuals(loess(dusk$hours[!dusk$filter]~
                                                         as.numeric(dusk$datetime[!dusk$filter]),span=0.1))),k=d)
  
  if(length(del.dawn)&gt;0) dawn$filter[!dawn$filter][del.dawn] &lt;- TRUE
  if(length(del.dusk)&gt;0) dusk$filter[!dusk$filter][del.dusk] &lt;- TRUE
  }
  
  if(plot){
  par(mfrow=c(2,1),mar=c(3,3,0.5,3),oma=c(2,2,0,0))
  plot(dawn$datetime[dawn$type==1],dawn$hours[dawn$type==1],pch="+",cex=0.6,xlab="",ylab="",yaxt="n")
  lines(dawn$datetime[!dawn$filter], predict(loess(dawn$hours[!dawn$filter]~as.numeric(dawn$datetime[!dawn$filter]),span=0.1)) , type="l")
  points(dawn$datetime[dawn$filter],dawn$hours[dawn$filter],col="red",pch="+",cex=1)
  axis(2,labels=F)
  mtext("Sunrise",4,line=1.2)
  
  plot(dusk$datetime[dusk$type==2],dusk$hours[dusk$type==2],pch="+",cex=0.6,xlab="",ylab="",yaxt="n")
  lines(dusk$datetime[!dusk$filter], predict(loess(dusk$hours[!dusk$filter]~as.numeric(dusk$datetime[!dusk$filter]),span=0.1)), type="l")
  points(dusk$datetime[dusk$filter],dusk$hours[dusk$filter],col="red",pch="+",cex=1)
  axis(2,labels=F)
  legend("bottomleft",c("Outside filter","Inside filter"),pch=c("+","+"),col=c("black","red"),
         bty="n",cex=0.8)
  mtext("Sunset",4,line=1.2)
  mtext("Time",1,outer=T)
  mtext("Sunrise/Sunset hours (rescaled)",2,outer=T)
  }
  all &lt;- rbind(subset(dusk,filter),subset(dawn,filter))
  
  filter &lt;- rep(FALSE,length(tFirst))
  filter[tFirst%in%all$datetime | tSecond%in%all$datetime] &lt;- TRUE
  </code></pre>
  
  <p># original code:
    #return(!filter)</p>
  
  <p># altered code to return outliersreturn(del.dusk)
    # replace with code below to print outlying points</p>
  
  <pre><code>return(c("delete dawn",del.dawn,"delete dusk",del.dusk))
  
  }
  </code></pre>
  
  <p>APPLY FUNCTION
      loessFilter(log2$tFirst, log2$tSecond, type=1, k=4, plot=TRUE)</p>
  
  <p>remove the values - need to remove both sunrise and sunset curves</p>
  
  <pre><code>log2b&lt;-log2[-c(77,78,124,125),]
  length(log2$tFirst)
  length(log2b$tFirst)
  </code></pre>
  
  <p>repeat function to see if the values have gone</p>
  
  <pre><code>loessFilter(log2b$tFirst, log2b$tSecond, type=1, k=4, plot=TRUE)
  </code></pre>
  
  <p>outliers still there!!</p>
  
  <p>HERE ARE THE DATA:
  <a href="http://www.4shared.com/file/jxVuTsVHce/002_geolight.html" rel="nofollow">http://www.4shared.com/file/jxVuTsVHce/002_geolight.html</a></p>
  
  <p>A bit too long to post the full data here and the example won't work with a dummy dataset :)</p>