FIREFEEDER REPORT
Firefeeder data analysis
First, let’s look into a per-animal summmary
setwd("C:/Users/Admin/OneDrive/Documents/phenomics/FIRE")
#setwd("C:/Users/Marti/OneDrive/Documents/phenomics/FIRE")
invisible(library(tidyverse))
invisible(library(DT))
invisible(library(lubridate))
invisible(library(lme4))
invisible(library(broom))
fire_file<-"FIRE_4.25.18.csv"
fire_data<-read_csv(fire_file)
#fire_data
key<-read_csv("FIRE_key.csv")
#key
#zerotags
filter(fire_data%>%arrange(-Consumed),Tag=="000000000000000")## # A tibble: 696 x 10
## Location Tag Date Entry Exit `Ent Wt` `Ext Wt` Consumed Weight
## <int> <chr> <chr> <tim> <tim> <dbl> <dbl> <dbl> <dbl>
## 1 1 00000000~ 4/18/~ 04:14 04:15 0.877 0.577 0.300 0.
## 2 1 00000000~ 4/20/~ 23:46 23:46 0.981 0.685 0.296 0.
## 3 2 00000000~ 4/8/2~ 23:18 23:18 0.887 0.620 0.267 0.
## 4 1 00000000~ 4/25/~ 01:08 01:08 0.941 0.687 0.254 0.
## 5 1 00000000~ 4/8/2~ 06:50 06:50 0.831 0.581 0.250 0.
## 6 1 00000000~ 4/24/~ 02:48 02:48 1.10 0.858 0.247 0.
## 7 1 00000000~ 4/23/~ 19:44 19:46 0.783 0.666 0.117 102.
## 8 1 00000000~ 4/16/~ 19:39 19:39 0.594 0.515 0.0790 0.
## 9 1 00000000~ 4/19/~ 03:37 03:37 0.578 0.500 0.0780 0.
## 10 1 00000000~ 4/21/~ 21:43 21:43 0.649 0.571 0.0780 0.
## # ... with 686 more rows, and 1 more variable: `Topup Amount` <dbl>fire_data<-filter(fire_data,Tag!="000000000000000")
ffd<-mutate(fire_data,visit_length=difftime(Exit,Entry,units = "secs"),ID=round(as.numeric(Tag)-1000*985152014209))
fdn<-filter(ffd,visit_length<=0)
ffd<- mutate(ffd,Date=mdy(Date))
ffd<-mutate(ffd,edt=ymd_hms(Date+lubridate::hms(Entry)))
ffd<-mutate(ffd,exdt=ymd_hms(Date+(Exit<Entry)+lubridate::hms(Exit)))
ffd<-mutate(ffd,visit_length=exdt-edt)
ffd<-mutate(ffd,cweight=ifelse(Weight==0,NA,Weight))
idx<-match(ffd$ID,key$RFID_tag)
ffd<-(mutate(ffd,earnotch=key$ear_notch[idx]))
ffd<-filter(ffd,Date>="2018-04-02")
print(min(ffd$edt))## [1] "2018-04-02 08:26:59 UTC"print(max(ffd$exdt))## [1] "2018-04-25 09:35:11 UTC"dim(ffd)## [1] 3911 16sum_per_animal<-ffd %>%group_by(ID)%>% summarize(Tag=Tag[1],med_weight=median(cweight,na.rm = T),
total_feed=sum(Consumed),
visits=sum(Consumed*0+1),
total_time=(seconds(sum(visit_length))),
early=min(edt),
time_to_first_visit=round(difftime(min(edt),min(ffd$edt),units="hours"),2))
sum_per_animal<-mutate(sum_per_animal,total_time=round(total_time/3600,2))
idr<-match(key$RFID_tag,sum_per_animal$ID)
animal_report<-cbind(key,sum_per_animal[idr,-1])
animal_report<-animal_report[order(animal_report$total_time),]
datatable(select(animal_report,unit,Tag,ear_notch,med_weight,total_feed,visits,total_time,time_to_first_visit),
rownames=F,filter="top",
options = list(pageLength = 25,columnDefs = list(list(className = 'dt-center', targets = 0:7))))Relationship between time in feeder and feed use.
ggplot(ffd,aes(x=visit_length,y=Consumed))+geom_point()+geom_smooth()
ggplot(ffd,aes(x=visit_length,y=Consumed,color=earnotch))+geom_point()+geom_smooth(method="lm",se=F)
feedr<-mutate(ffd,kmin=visit_length/3600)%>%group_by(earnotch)%>%do(fit=lm(Consumed~kmin,data=.))
datatable(bind_cols(feedr%>%tidy(fit)%>%filter(term=="kmin")%>%select(earnotch,estimate,std.error,p.value),select(feedr%>%glance(fit),r.squared))%>%select(earnotch,estimate,p.value,r.squared),
rownames=F,
options = list(pageLength = 25,columnDefs = list(list(className = 'dt-center', targets = 0:3))),
filter="top"
)%>%formatRound(columns = c("estimate","p.value","r.squared"))Weight gain and feed efficiency
peso<-mutate(ffd,tod=as.numeric((edt-min(edt))/(3600*24)))%>%select(Location, earnotch, tod,cweight)
ggplot(peso,aes(x=tod,y=cweight,color=earnotch))+geom_point()+geom_smooth(se=F,method="lm")+xlab("Days")+ylab("Weight in Kg")+ylim(0,90)## Warning: Removed 1327 rows containing non-finite values (stat_smooth).## Warning: Removed 1327 rows containing missing values (geom_point).## Warning: Removed 62 rows containing missing values (geom_smooth).
gp<-peso%>%group_by(earnotch)%>%filter(!is.na(earnotch))%>%do((lm(cweight~tod,data=.) %>%summary()%>%coef()%>%round(x =.,3) %>%
as_data_frame())[2,c(1,2,4)])
#datatable(gp)
colnames(gp)<-c("earnotch","ADG","SE","P-val")
animal_report<-mutate(animal_report,
ADFI=round(total_feed/as.numeric(max(ffd$exdt)-min(ffd$edt)),3))
idx<-match(gp$earnotch,animal_report$ear_notch)
wss<-bind_cols(gp,animal_report[idx,c("ADFI","visits")])%>%mutate(feed_eff=ADG/ADFI)%>%arrange(feed_eff)%>%mutate(feed_eff=round(feed_eff,3))
datatable(wss,filter="top",rownames=F,
options = list(pageLength = 25,columnDefs = list(list(className = 'dt-center', targets = 0:6))))#Kg Feed per Kg of Gain
print(1/mean(wss$feed_eff))## [1] 2.786162Observation: the the column feed efficiency is kg of ADWG per kg of ADFI. The average Kg feed per Kg of gain is more or less on target with the expected values. Nut there are two very low performing animals that bring it down.
Behavior of dyads of pigs
A question we can as is what happens with a a pig’s feeder visit when another pig follows it quickly after. Is the first visit cut short? To answer this question we can compute the lag between the current visit and the next visit and record the pig of the current visit and the pig of the following visit.
gpd<-group_by(ffd,Location)
fdf<-function(dataset){
x<-dataset$edt[-1]-dataset$exdt[-nrow(dataset)]
print(sum(x<0))
ds<-tibble(Location=dataset$Location[-1],lag=x,Date=dataset$Date[-1],
Entry=dataset$edt[-nrow(dataset)],Exit=dataset$exdt[-nrow(dataset)],an1=dataset$earnotch[-nrow(dataset)],
Consumed=dataset$Consumed[-nrow(dataset)],preInt=dataset$visit_length[-nrow(dataset)],
an2=dataset$earnotch[-1])
ds<-ds[x>=0,]
return(ds)
}
vv<-fdf(filter(ffd,!is.na(earnotch)))## [1] 1vv<- filter(vv,Entry>="2018-04-02 04:30:00 UTC")%>%mutate(tonext=as.numeric(lag),current=as.numeric(preInt))
datatable(vv%>%select(Location,Entry,an1,current,Consumed,an2,tonext),rownames=F,filter="top",
options = list(pageLength = 5,columnDefs = list(list(className = 'dt-center', targets = 0:6))))#group_by(vv,an1,an2)%>%summarize(intensity=mean(preInt),freq=length(preInt))%>%ggplot(aes(as.factor(freq),intensity))+geom_boxplot()+geom_jitter(width = 0.1)+xlab('number of interactions')+ylab("average lag time")+ggtitle("Average Lag between consecutive visits and number of visits for a dyad of pigs")Does the identity of the pig following another one shortens up the previous visit? To answer this, I divided all visits in two groups: those that happen immediately after the current visit and those that took longer. As criteria for “immediate” I used 60 seconds. This need much improvement.
not immediate visits (more than 4 minutes after previous one ended)
trh=240
print(trh) ## [1] 240mtr=trh
tr=Inf
anova(lm(current~an1+an2+an1*an2,data=filter(vv,(tonext<tr)&(tonext>mtr))))## Analysis of Variance Table
##
## Response: current
## Df Sum Sq Mean Sq F value Pr(>F)
## an1 23 49987889 2173386 6.9367 <2e-16 ***
## an2 22 9707194 441236 1.4083 0.1062
## an1:an2 182 39938296 219441 0.7004 0.9964
## Residuals 352 110287437 313317
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Consumed~an1+an2+an1*an2,data=filter(vv,(tonext<tr)&(tonext>mtr))))## Analysis of Variance Table
##
## Response: Consumed
## Df Sum Sq Mean Sq F value Pr(>F)
## an1 23 10.6566 0.46333 6.8838 < 2e-16 ***
## an2 22 2.3388 0.10631 1.5795 0.04836 *
## an1:an2 182 10.1036 0.05551 0.8248 0.92795
## Residuals 352 23.6921 0.06731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mf<-lmer(current~(1|an1)+(1|an2)+(1|an1:an2),data=filter(vv,(tonext<tr)&(tonext>mtr)))
m0<-lmer(current~(1|an1),data=filter(vv,(tonext<tr)&(tonext>mtr)))
anova(m0,mf)## Data: filter(vv, (tonext < tr) & (tonext > mtr))
## Models:
## m0: current ~ (1 | an1)
## mf: current ~ (1 | an1) + (1 | an2) + (1 | an1:an2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 8988.1 9001.2 -4491.0 8982.1
## mf 5 8987.7 9009.5 -4488.9 8977.7 4.3417 2 0.1141nm<-rownames(ranef(mf)$an2)
mm<-cbind(ranef(mf)$an1[nm,],ranef(mf)$an2)Immediate visits
mtr=0
tr=trh
anova(lm(current~an1+an2+an1*an2,data=filter(vv,(tonext<tr)&(tonext>mtr))))## Analysis of Variance Table
##
## Response: current
## Df Sum Sq Mean Sq F value Pr(>F)
## an1 23 153295480 6665021 28.2239 < 2.2e-16 ***
## an2 22 26696909 1213496 5.1387 6.272e-14 ***
## an1:an2 221 73054246 330562 1.3998 0.0001514 ***
## Residuals 3057 721905258 236148
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Consumed~an1+an2+an1*an2,data=filter(vv,(tonext<tr)&(tonext>mtr))))## Analysis of Variance Table
##
## Response: Consumed
## Df Sum Sq Mean Sq F value Pr(>F)
## an1 23 45.322 1.97050 30.9189 < 2.2e-16 ***
## an2 22 7.307 0.33214 5.2116 3.306e-14 ***
## an1:an2 221 20.541 0.09295 1.4584 2.306e-05 ***
## Residuals 3057 194.827 0.06373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mf<-lmer(current~(1|an1)+(1|an2)+(1|an1:an2),data=filter(vv,(tonext<tr)&(tonext>mtr)))
m0<-lmer(current~(1|an1),data=filter(vv,(tonext<tr)&(tonext>mtr)))
anova(m0,mf)## Data: filter(vv, (tonext < tr) & (tonext > mtr))
## Models:
## m0: current ~ (1 | an1)
## mf: current ~ (1 | an1) + (1 | an2) + (1 | an1:an2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 50819 50837 -25407 50813
## mf 5 50761 50792 -25376 50751 61.799 2 3.806e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mf1<-lmer(current~(1|an1)+(1|an2),data=filter(vv,(tonext<tr)&(tonext>mtr)))
anova(mf1,mf)## Data: filter(vv, (tonext < tr) & (tonext > mtr))
## Models:
## mf1: current ~ (1 | an1) + (1 | an2)
## mf: current ~ (1 | an1) + (1 | an2) + (1 | an1:an2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mf1 4 50762 50786 -25377 50754
## mf 5 50761 50792 -25376 50751 2.7355 1 0.09814 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mf1,m0)## Data: filter(vv, (tonext < tr) & (tonext > mtr))
## Models:
## m0: current ~ (1 | an1)
## mf1: current ~ (1 | an1) + (1 | an2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 50819 50837 -25407 50813
## mf1 4 50762 50786 -25377 50754 59.063 1 1.527e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#nm<-rownames(ranef(mf)$an2)
#mm<-cbind(ranef(mf)$an1[nm,],ranef(mf)$an2)
#plot(mm)As we can see, the identity of the current pig (an1) always influences the length of the visit (some pigs take longer at the feeder than others). But the identity of the following pig only is associated with the length of the current visit if it follows it quickly after. This does not mean that the following pig always shortens the visit, it may well be that pigs go to feed together sometimes and the effect of the follower could be just confounded with the effect of the current pig.
Feeder Use.
A worry I have is if we are having too much competition due to overloading the feeder with pigs. first I computed the total use time per day. Then I computed it by hour.
sum_per_feeder<-ffd%>%group_by(Location,Date)%>%summarize(totaluse=round(sum(visit_length)/3600,2),Feed=sum(Consumed))
datatable(sum_per_feeder,rownames = F,options = list(pageLength = 25,columnDefs = list(list(className = 'dt-center', targets = 0:3))))lg<-mutate(vv,lng=1+(tonext>=tr))
ct<-group_by(vv,hour(Exit))%>%summarize(lgs=3600-mean(tonext))
colnames(ct)<-c("hr","lgs")
per<-24
rg<-3600/2
reslm <- lm(ct$lgs ~ sin(2*pi/per*ct$hr)+cos(2*pi/per*ct$hr))
ft<-fitted(reslm)
ggplot(ct,aes(hr,lgs))+geom_point()+xlab("Hour of the day")+ylab("Feeder use per hour in seconds")
#+geom_smooth(span=.5,se=F)+geom_hline(yintercept=3600)+geom_line(aes(y=ft/max(ft)*3600,color="red"))+xlab("Hour of the day")+ylab("Feeder use per hour in seconds")We need to discuss the weight gain tables and the results of the feeder use.