FIREFEEDER REPORT 2
Firefeeder data analysis
First, let’s look into a per-animal summmary
setwd("C:/Users/Admin/OneDrive/Documents/phenomics/FIRE/trial2")
#setwd("C:/Users/Marti/OneDrive/Documents/phenomics/FIRE/trial2")
invisible(library(tidyverse))
invisible(library(DT))
invisible(library(lubridate))
invisible(library(lme4))
invisible(library(broom))
fire_file<-"FIRE_8.23.18.csv"
fire_data<-read_csv(fire_file)
#fire_data
key<-read_csv("key2.csv")
#key
fire_data<-filter(fire_data,Tag!="128034494483425")
#fire_data
#zerotags
zd<-filter(fire_data%>%arrange(-Consumed),Tag=="000000000000000")
zd<- mutate(zd,Date=mdy(Date))
zd<-mutate(zd,edt=ymd_hms(Date+lubridate::hms(Entry)))
zd<-mutate(zd,exdt=ymd_hms(Date+(Exit<Entry)+lubridate::hms(Exit)))
zd<-mutate(zd,visit_length=exdt-edt)
zd<-mutate(zd,cweight=ifelse(Weight==0,NA,Weight))
zd<-filter(zd,Date>="2018-06-13"&Date<"2018-07-10")
zd<-filter(zd,visit_length>5,Consumed>0.1)
#filter(zd,)
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-06-13")
ffd<-filter(ffd,!is.na(earnotch))
print(min(ffd$edt))## [1] "2018-06-13 00:15:14 UTC"print(max(ffd$exdt))## [1] "2018-08-21 14:17:21 UTC"dim(ffd)## [1] 18982 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)%>%filter(!is.na(earnotch)))%>%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 3162 rows containing non-finite values (stat_smooth).## Warning: Removed 3162 rows containing missing values (geom_point).## Warning: Removed 30 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] 1.588957Observation: 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] 2vv<- 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 31 54012388 1742335 12.0901 < 2.2e-16 ***
## an2 29 22786759 785750 5.4523 < 2.2e-16 ***
## an1:an2 290 86850975 299486 2.0781 < 2.2e-16 ***
## Residuals 1742 251044842 144113
## ---
## 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 31 3.431 0.110693 5.3823 < 2.2e-16 ***
## an2 29 1.700 0.058637 2.8511 7.205e-07 ***
## an1:an2 290 8.645 0.029812 1.4496 7.007e-06 ***
## Residuals 1742 35.826 0.020566
## ---
## 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 31308 31325 -15651 31302
## mf 5 31270 31298 -15630 31260 42.02 2 7.508e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nm<-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 31 121936327 3933430 20.7143 < 2.2e-16 ***
## an2 30 64104522 2136817 11.2530 < 2.2e-16 ***
## an1:an2 443 213577834 482117 2.5389 < 2.2e-16 ***
## Residuals 16377 3109818671 189889
## ---
## 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 31 15.70 0.50659 20.1970 < 2.2e-16 ***
## an2 30 8.48 0.28264 11.2683 < 2.2e-16 ***
## an1:an2 443 23.58 0.05324 2.1225 < 2.2e-16 ***
## Residuals 16377 410.78 0.02508
## ---
## 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 254163 254186 -127079 254157
## mf 5 253815 253854 -126903 253805 351.91 2 < 2.2e-16 ***
## ---
## 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 253955 253986 -126974 253947
## mf 5 253815 253854 -126903 253805 142.11 1 < 2.2e-16 ***
## ---
## 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 254163 254186 -127079 254157
## mf1 4 253955 253986 -126974 253947 209.8 1 < 2.2e-16 ***
## ---
## 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.