1.1. Mixed model on visit length at the feeder

\[y=X\beta + Zu + e\] Where, \(y\) is a \(n x 1\) vector of visit length at the feeder (minutes), \(X\beta\) are the fixed effects, \(location-trial \ (12)\), \(hour \ entry \ at \ the \ feeder \ (23)\) and \(animal \ median \ weigth\) as covariate; \(Z\) is a \(n\) x \(q\) desing matrix (\(q\) is the number of pigs) relates records in \(y\) to the random vector of additive genetic effects \(u\) (\(q\) x \(1\)); \(e\) (\(n\) x \(1\)) is the random residuals vector.

tn<-60
trials60.data<-trials6.data%>%filter(to_next<=tn)
dim(trials60.data)

## [1] 50561    22

m0.60s<-lmer(visit.length ~ Loc_trial -1 + hour_entry+ wt_median+
               (1|Eartag_trial),data=trials60.data)
a<-summary(m0.60s)
a$coefficients%>%as.data.frame()%>%kable(caption = "Estimates of fixed effects")%>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left", font_size = 12)

a0<-as.matrix(as.data.frame(VarCorr(m0.60s))[,4],1,2)
colnames(a0)<-"Variance"
rownames(a0)<-c("Eartag_M0","Residual_M0")

1.2. Mixed model: Location + Hours + median weight + Eartag + Follower

\[y=X\beta + Zu + Z_fa_f + e\] Where, \(y\) is a \(n x 1\) vector of visit length at the feeder (minutes), \(X\beta\) are the fixed effects, \(location-trial \ (12)\), \(hour \ entry \ at \ the \ feeder \ (23)\) and \(animal \ median \ weigth\) as covariate; \(Z\) is a \(n\) x \(q\) desing matrix (\(q\) is the number of pigs) relates records in \(y\) to the random vector of additive genetic effects \(u\) (\(q\) x \(1\)); \(Z_f\) is the design matrix of the next individual that visited the feeder, named \(followers\), relating to \(y\) with the random vector effects \(a_f\) (\(q\) x \(1\)) and \(e\) (\(n\) x \(1\)) is the random residuals vector.

m1.60s<-lmer(visit.length ~ Loc_trial -1 + hour_entry+ wt_median+
           (1|Eartag_trial)+(1|follower_trial),data=trials60.data)
b<-summary(m1.60s)
b$coefficients%>%as.data.frame()%>%kable(caption = "Estimates of fixed effects")%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, 
              position = "left",font_size = 12)

b0<-as.matrix(as.data.frame(VarCorr(m1.60s))[,4],1,3)
colnames(b0)<-"Variance"
rownames(b0)<-c("Eartag_M1","Follower_M1","Residual_M1")

1.3. Estimated variance components

rbind(a0,b0)%>%kable(caption = "Estimated of variance components")%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left",font_size = 12)

1.4. Model comparison

anova(m0.60s,m1.60s)

## refitting model(s) with ML (instead of REML)

## Data: trials60.data
## Models:
## m0.60s: visit.length ~ Loc_trial - 1 + hour_entry + wt_median + (1 | 
## m0.60s:     Eartag_trial)
## m1.60s: visit.length ~ Loc_trial - 1 + hour_entry + wt_median + (1 | 
## m1.60s:     Eartag_trial) + (1 | follower_trial)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m0.60s 38 332624 332960 -166274   332548                             
## m1.60s 39 331643 331988 -165783   331565 983.12      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5. Daily Feed Intake by Animal in trial day

The animal average daily intake was calculated as:

\[AFI_{i} = \frac{\sum_{j=1}^n X_j}{total\ days\ trial} \]

Where \(X_{j}\) is the \(jth\) consumed of an individual (\(j=1,..n\))

Feed.intake<-trials6.data%>%group_by(Loc_trial,Eartag_trial)%>%
                            summarize(animal.intake=sum(Consumed),
                               total.days=(max(trial.day)-min(trial.day)),
                               Average.daily.intake=animal.intake/total.days)%>%
          ungroup()%>%as.data.frame()

rownames(Feed.intake)<-Feed.intake$Eartag_trial
#----------
ggplot(Feed.intake, aes(x=Loc_trial, y=Average.daily.intake))+geom_boxplot()+
  ggtitle("Average Daily Feed Intake (kg) by Trial") + 
  xlab("Location by Trial")+ ylab("Average Daily Feed Intake (Kg)")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none") 

#----------

lmod.feed.intake<-lm(Average.daily.intake ~ Loc_trial-1, data = Feed.intake)
summary(lmod.feed.intake)$coefficients%>%as.data.frame()%>%kable(caption = "Estimated feed intake by location")%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left",font_size = 12)

resi.feed.intake<-as.data.frame(lmod.feed.intake$residuals)%>%
  mutate(., feed.intake.fitloc=lmod.feed.intake$residuals,
         Eartag_trial=rownames(.))%>%select(feed.intake.fitloc,Eartag_trial)
Feed.intake<-left_join(Feed.intake,resi.feed.intake, by="Eartag_trial")
Feed.intake$Loc_trial<-as.character(Feed.intake$Loc_trial)
# output file:weight gain, feed intake, Blups 
ADGs.FI.BLUPs<-left_join(blps.gr,Feed.intake, by= c("Eartag_trial","Loc_trial", "total.days"))

6. Correlation between blups and Weight gain, feed intake

\(EarBLUP = animal\ ramdom \ effect\ of \ visit \ length \ at \ the \ feeder\)

\(FollBLUP = animal\ ramdom \ effect\ as \ next \ individual \ that \ visited \ the \ feeder\)

\(ADG= ADG_i\)

\(ADG.locmean=ADG_i-\hat{\ loc:trial}\)

\(weight.gain.qr= estimated\ weight\ gain\ by\ median\ regression\ model\)

\(trial.day.rr = trial\ day\ random\ effect\ of\ the\ ith\ animal\)

\(Weight.gain.rr = group\ slope + trial.day.rr\)

\(Average.daily.intake= AFI_i\)

\(feed.intake.fitloc = AFI_i- \hat{\ loc:trial}\)

ADGs.FI.BLUPs%>%select(EarBLUP,FollBLUP,ADG,ADG.locmean,
                       weight.gain.qr, trial.day.rr, Weight.gain.rr,
                       Average.daily.intake, feed.intake.fitloc)%>%cor()%>%
  kable()%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F,
                          position = "left",font_size = 12)

7. Visit length by animal

visitID.mean<-trials6.data%>%group_by(Eartag_trial, Loc_trial)%>%
  summarize(vl=sum(visit.length))%>%ungroup()%>%as.data.frame()
rownames(visitID.mean)<-visitID.mean$Eartag_trial
summary(visitID.mean$vl.sum)

## Length  Class   Mode 
##      0   NULL   NULL

vlm1<-lm(vl ~ Loc_trial-1, data = visitID.mean)
summary(vlm1)$coefficients%>%as.data.frame()%>%kable(caption = "Mean estimated by location-trial group")%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left",font_size = 12)

visitID.mean$Loc_trial<-as.character(visitID.mean$Loc_trial)

7.2. Random regression model on Visit lenght

The random regression model on visit length in each trial in matrix notation is:

\[Y_{j}= X_j\beta + Z_ju_{j} + \ e_{j}\]

Where:

\(Y_{ij}\) is a matrix \(n_j\)x\(1\) with all observation of visit length total for each day, in \(jth\) trial

\(X_j\) is the design matrix \(n_j\)x\(2\) of fixed effects for each group (location-trial, trial day)

\(\beta\) is the fixed effects vector, with \(\beta=(\beta_{0}, \beta_{1})'\)

\(Z_j\) is the design matrix \(n_j\)x\(2\) of random effects for each group (Animal, trial day)

\(u_j\) is the random effects vector, with \(u=(u_{0j}, u_{1j})'\)

\(e_{ij}\) is the error vector

and the the \(ith\) row (\(i = 1,...n_j\)) of the response vector is:

\[y_{ij}= \beta_0 + \beta_1x_{ij} + u_0 + u_{1j}x_{ij} + e_{ij}\]

The visit length by animal is

\(Visit.length_{pig}= group\ slope + trial.day.rr = (\beta_1 + u_{1j})x_{ij}\)

vlmday<-trials6.data%>%group_by(Loc_trial,Eartag_trial, trial.day)%>%
  summarize(visit.mean.day=sum(visit.length))%>%ungroup()
vlm.day.lm1<-lmer(visit.mean.day ~ 1 + Loc_trial*trial.day +
               (trial.day| Eartag_trial), data = vlmday,
             control=lmerControl(optCtrl=list(maxfun=50000)))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00497934
## (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

# Estimates
b<-summary(vlm.day.lm1)
b$coefficients%>%as.data.frame()%>%kable(caption = "Estimates of fixed effects")%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, 
              position = "left",font_size = 12)

# Variances
print(VarCorr(vlm.day.lm1), comp="Variance")

##  Groups       Name        Variance  Corr  
##  Eartag_trial (Intercept) 497.02632       
##               trial.day     0.10954 -0.591
##  Residual                 359.34987

# fixed effects
fxeff<-as.matrix(summary(vlm.day.lm1)$coef[,1])
fxeff<-tibble(Loc_trial=rownames(fxeff)[1:12],
              intercept.vl=fxeff[1:12,],slope.vl=fxeff[13:24,])%>%
        mutate(.,Loc_trial=str_sub(Loc_trial,start = 10,end = 21))
fxeff[1,1]<-"Loc_1_t_1"
# intercetpt group
fxeff$intercept.vl[-1]<-fxeff$intercept.vl[-1]+fxeff$intercept.vl[1] # sum means
# slope group
fxeff$slope.vl[-1]<-fxeff$slope.vl[-1]+fxeff$slope.vl[1] 
fxeff%>%head()%>%kable()%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left", font_size = 12)

# Blups visit length rr
blps.vl<-ranef(vlm.day.lm1)
blps.vl<-blps.vl[[1]]%>%mutate(ID=rownames(blps.vl[[1]]))%>%
  select(ID,  "(Intercept)", "trial.day")
colnames(blps.vl)<-c("ID", "EarBlup.vl.rr","vl.day")
blps.vl<-left_join(ADGs.FI.BLUPs,blps.vl, by = c("Eartag_trial"="ID"))
blps.vl$Loc_trial<-as.character(blps.vl$Loc_trial)
blps.vl<-left_join(blps.vl,fxeff,"Loc_trial")
output.vl<-mutate(blps.vl,Visit.length=vl.day+slope.vl)

7.3. Correlation between blups, visit length, weight gain, feed intake

\(EarBLUP = animal\ ramdom \ effect\ of \ visit \ length \ at \ the \ feeder\)

\(FollBLUP = animal\ ramdom \ effect\ as \ next \ individual \ that \ visited \ the \ feeder\)

\(trial.day.rr = trial\ day\ random\ effect\ of\ the\ ith\ animal\ on\ weight\ gain\)

\(feed.intake.fitloc = AFI_i- \hat{\ loc:trial}\)

\(visitl.fitloc = Vl_i-\hat{\ loc:trial}\)

\(vl.day = trial\ day\ random\ effect\ of\ the\ ith\ animal\ on\ visit\ length\ at\ the\ feeder\)

\(Visit.length_{pig}= group\ slope + vl.day = (\beta_1 + u_{1j})x_{ij}\)

kable(e1, caption = " Estimates correlation between BLUP's,visit length, weight gain and feed intake ") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, 
                position = "left",font_size = 14) %>%
  add_header_above(c(" " = 1, "Eartag Blup" = 2, "Follower Blup" = 2))

#output.vl%>%select("EarBLUP","FollBLUP","trial.day.rr","feed.intake.fitloc",
 #                  "visitl.fitloc", "Visit.length","vl.day")%>%cor()