1. Mixed model on visit length at the feeder

1.1. Models with set 1: next visit <= 60 sc 58255 records

1.1.1. Mixed model: Location + Hours + median weight + Eartag

\[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 \ (14)\), \(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.

dim(trials.data)
## [1] 74415    22
tn<-60
trials60.data<-trials.data%>%filter(to_next<=tn)
dim(trials60.data)
## [1] 58255    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)
Estimates of fixed effects
Estimate Std. Error df t value Pr(>|t|)
Loc_trialLoc_1_t_1 13.3523943 0.9367143 145.5383 14.2545007 0.0000000
Loc_trialLoc_1_t_2 7.1239060 0.9547475 132.0047 7.4615604 0.0000000
Loc_trialLoc_1_t_3 10.1277485 0.9536971 131.4240 10.6194597 0.0000000
Loc_trialLoc_1_t_4 8.1740262 1.1128355 128.8692 7.3452240 0.0000000
Loc_trialLoc_1_t_5 9.0033900 1.1182379 131.3855 8.0514083 0.0000000
Loc_trialLoc_1_t_6 4.8740442 1.0477011 128.1404 4.6521323 0.0000081
Loc_trialLoc_1_t_7 10.9622969 1.0517354 130.1224 10.4230561 0.0000000
Loc_trialLoc_2_t_1 10.3959227 0.9316799 142.4271 11.1582559 0.0000000
Loc_trialLoc_2_t_2 7.9545331 0.9556811 132.5214 8.3234177 0.0000000
Loc_trialLoc_2_t_3 14.5941646 0.9553291 132.3243 15.2765837 0.0000000
Loc_trialLoc_2_t_4 6.8140992 1.1108549 127.9559 6.1341037 0.0000000
Loc_trialLoc_2_t_5 9.7234563 1.1204119 132.4073 8.6784655 0.0000000
Loc_trialLoc_2_t_6 6.8090292 1.0485406 128.5513 6.4938153 0.0000000
Loc_trialLoc_2_t_7 8.3649306 1.1132025 129.0404 7.5142940 0.0000000
hour_entry1 0.2247501 0.2811512 58098.3547 0.7993922 0.4240663
hour_entry2 0.0048850 0.2820536 58098.3767 0.0173196 0.9861817
hour_entry3 -0.1438705 0.2836057 58100.1968 -0.5072907 0.6119528
hour_entry4 -0.3333335 0.2744610 58100.3324 -1.2145022 0.2245609
hour_entry5 -0.5985261 0.2505983 58102.1352 -2.3883886 0.0169256
hour_entry6 -0.5939477 0.2301661 58102.3352 -2.5805181 0.0098676
hour_entry7 -1.1285361 0.2179817 58102.8096 -5.1772057 0.0000002
hour_entry8 -2.3064472 0.2109444 58101.9572 -10.9339100 0.0000000
hour_entry9 -1.0038285 0.2157097 58102.6313 -4.6536086 0.0000033
hour_entry10 -0.1038588 0.2192123 58101.8554 -0.4737815 0.6356575
hour_entry11 0.3988215 0.2194233 58101.6679 1.8175898 0.0691320
hour_entry12 0.6913416 0.2179567 58101.9668 3.1719218 0.0015151
hour_entry13 0.4050796 0.2157590 58102.3490 1.8774636 0.0604596
hour_entry14 0.5516624 0.2165836 58101.5035 2.5471102 0.0108645
hour_entry15 0.8032194 0.2187920 58101.7974 3.6711546 0.0002417
hour_entry16 1.3698887 0.2272425 58104.6077 6.0283122 0.0000000
hour_entry17 1.7253099 0.2431562 58103.1097 7.0954800 0.0000000
hour_entry18 1.2274118 0.2489077 58101.4001 4.9311923 0.0000008
hour_entry19 0.7756071 0.2518216 58100.6252 3.0799867 0.0020711
hour_entry20 0.4299856 0.2498648 58101.5188 1.7208730 0.0852792
hour_entry21 0.2550130 0.2550467 58099.3403 0.9998678 0.3173786
hour_entry22 -0.0986655 0.2535773 58098.8394 -0.3890944 0.6972077
hour_entry23 0.3999527 0.2634583 58099.0008 1.5180874 0.1289978
wt_median 0.0125074 0.0016885 58104.6209 7.4073793 0.0000000 
a0<-as.matrix(as.data.frame(VarCorr(m0.60s))[,4],1,2)
colnames(a0)<-"Estimate"
rownames(a0)<-c("Eartag M1","Residual M1")

1.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)
Estimates of fixed effects
Estimate Std. Error df t value Pr(>|t|)
Loc_trialLoc_1_t_1 13.2120009 0.9878285 179.0739 13.3747924 0.0000000
Loc_trialLoc_1_t_2 6.9287284 1.0092550 164.0063 6.8651909 0.0000000
Loc_trialLoc_1_t_3 9.8877423 1.0084918 163.5081 9.8044846 0.0000000
Loc_trialLoc_1_t_4 8.0393686 1.1775105 160.7301 6.8274286 0.0000000
Loc_trialLoc_1_t_5 9.0030599 1.1836679 164.0954 7.6060689 0.0000000
Loc_trialLoc_1_t_6 4.8184240 1.1085663 159.8106 4.3465366 0.0000245
Loc_trialLoc_1_t_7 11.1022864 1.1129202 162.3259 9.9758152 0.0000000
Loc_trialLoc_2_t_1 10.6354044 0.9835247 175.9664 10.8135608 0.0000000
Loc_trialLoc_2_t_2 7.8608718 1.0101263 164.5727 7.7820684 0.0000000
Loc_trialLoc_2_t_3 14.5895874 1.0100358 164.5092 14.4446243 0.0000000
Loc_trialLoc_2_t_4 6.5528795 1.1754573 159.6173 5.5747489 0.0000001
Loc_trialLoc_2_t_5 9.4895420 1.1851368 164.9076 8.0071281 0.0000000
Loc_trialLoc_2_t_6 7.0019242 1.1095585 160.3818 6.3105501 0.0000000
Loc_trialLoc_2_t_7 8.4109268 1.1778816 160.9343 7.1407237 0.0000000
hour_entry1 0.2027794 0.2775645 57999.9377 0.7305669 0.4650467
hour_entry2 -0.0125857 0.2785423 58003.0862 -0.0451842 0.9639607
hour_entry3 -0.1604122 0.2801533 58009.1542 -0.5725871 0.5669265
hour_entry4 -0.4357640 0.2713935 58032.2158 -1.6056535 0.1083555
hour_entry5 -0.6268503 0.2481418 58051.7811 -2.5261773 0.0115338
hour_entry6 -0.5312345 0.2279737 58047.7314 -2.3302447 0.0197966
hour_entry7 -1.0093495 0.2157921 58043.0058 -4.6774176 0.0000029
hour_entry8 -2.1636716 0.2086747 58032.3930 -10.3686325 0.0000000
hour_entry9 -0.9194867 0.2134291 58038.3596 -4.3081601 0.0000165
hour_entry10 0.0192920 0.2167952 58034.6814 0.0889872 0.9290924
hour_entry11 0.5087501 0.2170289 58035.0124 2.3441581 0.0190734
hour_entry12 0.8505692 0.2155970 58032.9054 3.9451818 0.0000798
hour_entry13 0.5828707 0.2134000 58031.7878 2.7313526 0.0063094
hour_entry14 0.7219648 0.2141570 58027.7187 3.3711939 0.0007489
hour_entry15 0.8825843 0.2163289 58026.6482 4.0798259 0.0000451
hour_entry16 1.4031450 0.2248506 58038.9322 6.2403438 0.0000000
hour_entry17 1.7442151 0.2405581 58032.0740 7.2507024 0.0000000
hour_entry18 1.2473646 0.2461968 58025.5817 5.0665356 0.0000004
hour_entry19 0.7659805 0.2490722 58022.8211 3.0753349 0.0021036
hour_entry20 0.5147546 0.2469738 58019.2820 2.0842480 0.0371420
hour_entry21 0.3195946 0.2519284 58006.7265 1.2685930 0.2045914
hour_entry22 -0.0294404 0.2504649 58007.8570 -0.1175428 0.9064303
hour_entry23 0.4419855 0.2601454 58004.0054 1.6989942 0.0893256
wt_median 0.0133051 0.0016754 58063.9775 7.9414191 0.0000000 
b0<-as.matrix(as.data.frame(VarCorr(m1.60s))[,4],1,3)
colnames(b0)<-"Estimate"
rownames(b0)<-c("Eartag M2","Follower M2","Residual M2")

1.2. Models with set 1: next visit >= 600 sc 6258 records

1.2.1. Mixed model: Location + Hours + median weight + Eartag

tn<-600
trials600.data<-trials.data%>%filter(to_next>=tn)
dim(trials600.data)
## [1] 6258   22
m0.600s<-lmer(visit.length ~ Loc_trial -1 + hour_entry+ wt_median+
               (1|Eartag_trial),data=trials600.data)
a1<-as.matrix(as.data.frame(VarCorr(m0.600s))[,4],1,2)
colnames(a1)<-"Estimate"
rownames(a1)<-c("Eartag M1","Residual M1")

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

m1.600s<-lmer(visit.length ~ Loc_trial -1 + hour_entry+ wt_median+
               (1|Eartag_trial)+(1|follower_trial),data=trials600.data)
## boundary (singular) fit: see ?isSingular
b1<-as.matrix(as.data.frame(VarCorr(m1.600s))[,4],1,3)
colnames(b1)<-"Estimate"
rownames(b1)<-c("Eartag M2","Follower M2","Residual M2")

1.3. Estimated variance components

as.data.frame(cbind(rbind(a0,b0), rbind(a1,b1)))%>%kable(caption = "Estimated of variance components")%>%kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 12)%>%add_header_above(c("Variance "=1,"Subset <= 60 s"= 1, "Subset >= 600 s"=1))%>%
  row_spec(3:5, bold = T)
Estimated of variance components
Variance
Subset <= 60 s
Subset >= 600 s
Estimate Estimate
Eartag M1 9.419249 10.63967
Residual M1 43.048618 47.76711
Eartag M2 9.346853 10.63967
Follower M2 1.261856 0.00000
Residual M2 41.894731 47.76711

1.4. Model comparison

Subset <= 60 s

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 40 385127 385486 -192524   385047                             
## m1.60s 41 383848 384216 -191883   383766 1281.4      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Subset >= 600 s

anova(m0.600s,m1.600s)
## refitting model(s) with ML (instead of REML)
## Data: trials600.data
## Models:
## m0.600s: visit.length ~ Loc_trial - 1 + hour_entry + wt_median + (1 | 
## m0.600s:     Eartag_trial)
## m1.600s: visit.length ~ Loc_trial - 1 + hour_entry + wt_median + (1 | 
## m1.600s:     Eartag_trial) + (1 | follower_trial)
##         Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0.600s 40 42303 42573 -21112    42223                        
## m1.600s 41 42305 42582 -21112    42223     0      1          1

1.5. Plot Blup Eartag vs Blup Follower from estimated Subset <= 60 s

ggplot(BlupEF, aes(x=EarBLUP, y=FollBLUP))+
  geom_point(aes(color = EarBLUP, size = FollBLUP)) +
  geom_smooth(method = "lm",se = FALSE, color = "red")+
  scale_color_continuous(low = "black", high = "blue") +
  guides(color = FALSE, size = FALSE) +
  geom_point(aes(y = FollBLUP), shape = 1) +
  labs(title = "Direct BV vs Follower BV",
       subtitle = paste("r =",round(cor(BlupEF$EarBLUP,BlupEF$FollBLUP), 5), ";","p-val =",
                        round(cor.test(BlupEF$EarBLUP,BlupEF$FollBLUP)$p.val,5)))+
  theme(plot.subtitle = element_text(size = 8,face = "bold"))+
  xlab("Direct Breeding value")+ylab("Follower Breeding value")

2. Average Daily Gain (ADG) by Animal corected by the group mean

The Average daily gain \(ADG\) for each individual was calculated as:

\[ADG_i=\frac{\ median(Wt_{initial}) \ - median(Wt_{end})}{total\ trial \ days}\] Where:

\(ADG_i\) is the Average daily gain in kilograms of an individual \(i\) in any trial, with \(i=1,..,135\)

\(median(Wt_{initial})\) is the median weight in the start day of any trial.

\(median(Wt_{end})\) is the median weight in the end day of any trial.

\(total\ trial \ days\) are the number of total days of each trial.

ggplot(ADG, aes(x=Loc_trial, y=ADG))+geom_boxplot()+
  ggtitle("Average Daily Gain (kg) by Trial") + 
  xlab("Location by Trial")+ ylab("Average Daily Gain (Kg)")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none")

2.1. Linear model on Average daily Gain (ADG)

\[y=X\beta + e\] Where:

\(y\) is a \(n x 1\) vector of Average daily gain (\(ADG_i\) kg) of an individual

\(X\beta\) are the fixed effects, \(location-trial \ (14)\)

\(e\) (\(n\) x \(1\)) is the random residuals vector.

and \(y_i-x_i\hat{\beta}=\hat{e_i}\) is the ADG of the animal \(i\)

\(Average\ daily\ gain\ as\ deviation\ of\ group\)

# Average daily gain as deviation of group 
adg.lm<-lm(ADG ~ Loc_trial-1, data = ADG)
summary(adg.lm)$coefficients%>%as.data.frame()%>%kable()%>%kable_styling()
Estimate Std. Error t value Pr(>|t|)
Loc_trialLoc_1_t_1 0.8726667 0.0293656 29.71733 0
Loc_trialLoc_1_t_2 0.8155569 0.0306713 26.59019 0
Loc_trialLoc_1_t_3 1.1307487 0.0306713 36.86661 0
Loc_trialLoc_1_t_4 0.8993560 0.0359653 25.00618 0
Loc_trialLoc_1_t_5 1.0890625 0.0359653 30.28089 0
Loc_trialLoc_1_t_6 0.8795238 0.0339084 25.93819 0
Loc_trialLoc_1_t_7 1.3302254 0.0339084 39.22991 0
Loc_trialLoc_2_t_1 0.9107778 0.0293656 31.01515 0
Loc_trialLoc_2_t_2 0.8055718 0.0306713 26.26464 0
Loc_trialLoc_2_t_3 1.0510695 0.0306713 34.26878 0
Loc_trialLoc_2_t_4 1.0985078 0.0359653 30.54351 0
Loc_trialLoc_2_t_5 1.0332031 0.0359653 28.72775 0
Loc_trialLoc_2_t_6 0.9412698 0.0339084 27.75915 0
Loc_trialLoc_2_t_7 1.1788043 0.0359653 32.77612 0 
resi.adg<-as.data.frame(adg.lm$residuals)%>%
  mutate(., ADG.locmean=adg.lm$residuals,
         Eartag_trial=rownames(.))%>%select(ADG.locmean,Eartag_trial)
ADG<-left_join(ADG,resi.adg)
## Joining, by = "Eartag_trial"

3. Regression model for quantile 0.5 of the Weight by individual

The regression model for quantile level \(\tau\) of the Weight by individual is:

\[Q_{\tau}(y_i)= \beta_{0}(\tau) + \beta_1(\tau)x_{i1}\] \(\tau=0.5\) median quantile

\((y_i)\) is the \(i-th\) observation of the Weight of an animal

\(\beta_{0}\) is the intercept parameter

\(\beta_1\) is the slope

\(x_{i1}\) is the trial day

The goodness of fit for quantile regression by individual, which is analog to \(R^2\), was calculated following to Koenker and Machado (1999), and then

\[R^1(\tau)=1-\frac{\hat{V}(\tau)}{\tilde{V}(\tau)}\] Where:

\(R^1(\tau)\) is a local measure of goodness of fit for a particular quantile

\(\hat{V}(\tau)\) is a model with the intercept parameter only

\(\tilde{V}(\tau)\) is the unrestricted model

The Weight gain (\(weight.gain.qr\))

suppressMessages(library(quantreg))

# 3.1. Quantile regression function, to estimate ADG by trial day
qrf<-function(x){
  # 1. Regression 0.5 quantile  model,  full model
      rst1<-trials.data%>% filter(Eartag_trial==x)%>%rq(Weight~trial.day,data=.,tau=0.5)
       ADG<-rst1$coefficients[2]
       rho1<-rst1$rho
       
  # 2. Regression 0.5 quantile  model,  restricted model (intercept only)
       rst0<-trials60.data%>% filter(Eartag_trial==x)%>%rq(Weight~1,data=.,tau=0.5)
       rho0<-rst0$rho
       
  # 3. Estimate R^2 by individual
       R_sq.rq<- 1 - rho1/rho0
       
  # 4. Output
       out1<-tibble(ADG=ADG,R_sq.rq=R_sq.rq)
       
      return(out1)
   }  

# 3.2. Mean Average Daily Gain  Quantile Regression
mns<-unique(trials.data$Eartag_trial)
wg<-sapply(mns,qrf)
wg<-as.matrix(t(wg))
Wg<-tibble(id=rownames(wg)%>%map_df(as_tibble), wg.qr=wg[,1]%>%map_df(as_tibble),
           R.sq=wg[,2]%>%map_df(as_tibble))%>%
  mutate(id=id$value, weight.gain.qr=wg.qr$value,
         R.sq=R.sq$value)%>%select(id, weight.gain.qr, R.sq)

Wg%>%summary()%>%kable()%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left",font_size = 12)
id weight.gain.qr R.sq
Length:135 Min. :0.6000 Min. :0.2593
Class :character 1st Qu.:0.8724 1st Qu.:0.7705
Mode :character Median :1.0000 Median :0.8559
NA Mean :0.9948 Mean :0.8182
NA 3rd Qu.:1.1174 3rd Qu.:0.9014
NA Max. :1.6933 Max. :0.9338 
hist(Wg$weight.gain.qr, main = "Histogram Daily Weight Gain (kg)",
     xlab = "Weight gain (kg)", col = "seagreen")

hist(Wg$R.sq, main = bquote(R^2~ "of regression 0.5 quantile on Animal Weight"),
     xlab = bquote(R^2), col = "royalblue4")

# Animals with R square smaller and greater
Wg[Wg$R.sq<0.5,]
## # A tibble: 2 x 3
##   id       weight.gain.qr  R.sq
##   <chr>             <dbl> <dbl>
## 1 1361_t_2          0.812 0.259
## 2 1132_t_3          0.6   0.389
Wg[Wg$R.sq>0.925,]
## # A tibble: 3 x 3
##   id       weight.gain.qr  R.sq
##   <chr>             <dbl> <dbl>
## 1 1672_t_3           1.11 0.928
## 2 1413_t_3           1.12 0.934
## 3 1467_t_7           1.26 0.934
trials.data%>%filter(Eartag_trial%in%c("1361_t_2","1132_t_3","1672_t_3","1413_t_3","1467_t_7"))%>%
  ggplot(aes(x=trial.day, y=Weight)) +  geom_point()+ facet_wrap(.~Eartag_trial) + 
  geom_quantile(quantiles = .5, size=0.75, color="chartreuse") + 
  theme_light()+ 
  ggtitle("Body weight (kg) by Trial Day") + labs(subtitle = bquote(R^2~"= 0.259 1361_t_2 |" ~ R^2~"= 0.388 1132_t_3 |" ~ R^2~"= 0.927 1672_t_3 |" ~ R^2~"= 0.933 1413_t_3 |" ~R^2~"= 0.934 1467_t_7"))+theme(plot.subtitle = element_text(size = 8))+ xlab("Trial day")+ ylab("Weight (Kg)")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none") 
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x

ADG<-left_join(ADG, Wg, by=c("Eartag_trial"="id"))

outADG.Blups<-left_join(ADG,BlupEF, by=c("Eartag_trial"="ID"))

4. Random regression model on median weight

The random regression model on median weight 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 median weight for each day, in \(jth\) trial

\(X_j\) is the design matrix \(n_j\)x\(2\) of fixed effects for each group

\(\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

\(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}\]

# 4.1. Animal Median weight by day data
AmWT.day<-trials.data%>%group_by(Loc_trial,Eartag_trial,trial.day)%>%
  summarize(weightmd.day=median(Weight))%>%
  ungroup()%>%as.data.frame()

# 4.2. Ramdom regresion model of median weight 
mod_gr<-lmer(weightmd.day ~ 1+Loc_trial*trial.day+
       (trial.day| Eartag_trial), data = AmWT.day,
     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.00644689 (tol = 0.002, component 1)
# Estimates
b<-summary(mod_gr)
b$coefficients%>%as.data.frame()%>%kable(caption = "Estimates of fixed effects random regression model")%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left",font_size = 12)
Estimates of fixed effects random regression model
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 72.5398063 1.1991919 132.0501 60.4905723 0.0000000
Loc_trialLoc_1_t_2 -37.1575194 1.7072849 124.0953 -21.7640996 0.0000000
Loc_trialLoc_1_t_3 -47.7619757 1.7069756 124.0040 -27.9804680 0.0000000
Loc_trialLoc_1_t_4 -32.2842224 1.8656189 123.7695 -17.3048323 0.0000000
Loc_trialLoc_1_t_5 -28.6484515 1.8815158 128.0401 -15.2262616 0.0000000
Loc_trialLoc_1_t_6 -44.8311845 1.7991747 122.8974 -24.9176385 0.0000000
Loc_trialLoc_1_t_7 -41.9809875 1.8065985 124.9376 -23.2375857 0.0000000
Loc_trialLoc_2_t_1 -0.9829362 1.6960995 132.1078 -0.5795274 0.5632203
Loc_trialLoc_2_t_2 -35.5552698 1.7072654 124.0896 -20.8258603 0.0000000
Loc_trialLoc_2_t_3 -50.1831078 1.7061330 123.7608 -29.4133616 0.0000000
Loc_trialLoc_2_t_4 -44.9280852 1.8654167 123.7159 -24.0847444 0.0000000
Loc_trialLoc_2_t_5 -30.5662212 1.8815158 128.0401 -16.2455300 0.0000000
Loc_trialLoc_2_t_6 -49.6295499 1.7991191 122.8822 -27.5854716 0.0000000
Loc_trialLoc_2_t_7 -34.0101650 1.8686289 124.5699 -18.2005987 0.0000000
trial.day 0.9210570 0.0349499 166.4911 26.3536034 0.0000000
Loc_trialLoc_1_t_2:trial.day -0.0985481 0.0480311 135.9138 -2.0517551 0.0421135
Loc_trialLoc_1_t_3:trial.day 0.2009205 0.0480224 135.8112 4.1838884 0.0000511
Loc_trialLoc_1_t_4:trial.day -0.0059824 0.0522569 133.2139 -0.1144801 0.9090296
Loc_trialLoc_1_t_5:trial.day 0.1665349 0.0535436 146.8013 3.1102672 0.0022461
Loc_trialLoc_1_t_6:trial.day -0.0029619 0.0504111 132.4372 -0.0587552 0.9532356
Loc_trialLoc_1_t_7:trial.day 0.3247141 0.0508000 136.5659 6.3920054 0.0000000
Loc_trialLoc_2_t_1:trial.day -0.0253662 0.0494516 166.8194 -0.5129492 0.6086654
Loc_trialLoc_2_t_2:trial.day -0.1113809 0.0480302 135.9040 -2.3189737 0.0218880
Loc_trialLoc_2_t_3:trial.day 0.1299458 0.0479838 135.3801 2.7081166 0.0076402
Loc_trialLoc_2_t_4:trial.day 0.1689095 0.0522456 133.0994 3.2329870 0.0015446
Loc_trialLoc_2_t_5:trial.day 0.1377457 0.0535436 146.8013 2.5725883 0.0110863
Loc_trialLoc_2_t_6:trial.day 0.0252125 0.0504106 132.4321 0.5001423 0.6178049
Loc_trialLoc_2_t_7:trial.day 0.2776523 0.0524412 135.1008 5.2945435 0.0000005 
# Variances
print(VarCorr(mod_gr), comp="Variance")
##  Groups       Name        Variance  Corr 
##  Eartag_trial (Intercept) 15.863553      
##               trial.day    0.011734 0.229
##  Residual                  4.276034
# fixed effects intercept and slope by trial
fxeff<-as.matrix(summary(mod_gr)$coef[,1])
fxeff<-tibble(Loc_trial=rownames(fxeff)[1:14],
              intercept.rr=fxeff[1:14,],slope.rr=fxeff[15:28,])%>%
        mutate(.,Loc_trial=str_sub(Loc_trial,start = 10,end = 21))
fxeff[1,1]<-"Loc_1_t_1"
# intercept group
fxeff$intercept.rr[-1]<-fxeff$intercept.rr[-1]+fxeff$intercept.rr[1] # sum means
# slope group
fxeff$slope.rr[-1]<-fxeff$slope.rr[-1]+fxeff$slope.rr[1] 
fxeff%>%kable(caption = "Intercept and Slope by group")%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left",font_size = 12)
Intercept and Slope by group
Loc_trial intercept.rr slope.rr
Loc_1_t_1 72.53981 0.9210570
Loc_1_t_2 35.38229 0.8225089
Loc_1_t_3 24.77783 1.1219775
Loc_1_t_4 40.25558 0.9150746
Loc_1_t_5 43.89135 1.0875919
Loc_1_t_6 27.70862 0.9180951
Loc_1_t_7 30.55882 1.2457711
Loc_2_t_1 71.55687 0.8956908
Loc_2_t_2 36.98454 0.8096761
Loc_2_t_3 22.35670 1.0510028
Loc_2_t_4 27.61172 1.0899665
Loc_2_t_5 41.97359 1.0588026
Loc_2_t_6 22.91026 0.9462695
Loc_2_t_7 38.52964 1.1987092 
# Blups rr
blps.gr<-ranef(mod_gr)
blps.gr<-blps.gr[[1]]%>%mutate(ID=rownames(blps.gr[[1]]))%>%
  select(ID,  "(Intercept)", "trial.day"  )
colnames(blps.gr)<-c("ID", "Eartag_blup.rr","trial.day.rr")
blps.gr<-left_join(outADG.Blups,blps.gr, by = c("Eartag_trial"="ID"))
blps.gr$Loc_trial<-as.character(blps.gr$Loc_trial)
blps.gr<-left_join(blps.gr,fxeff,"Loc_trial")
blps.gr<-mutate(blps.gr,Weight.gain.rr=trial.day.rr+slope.rr)

4.1. Plot individual Weight gain estimated random regression and ADG_i model

$weight.gain.rr_{pig}= group slope + trial.day.rr = (1 + u{1j})x_{ij} $

ggplot(blps.gr, aes(x=ADG, y=Weight.gain.rr))+geom_point() + geom_abline()+
  labs(title = bquote( "Weight gain estimated by random regression and" ~ "ADG"[i] ~ "models" ) )+ xlab(bquote("ADG"[i]))+ ylab("Weight gain random regression")

5. Daily Feed Intake by Animal in trial day

The animal average daily feed 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<-trials.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") 

#----------

5.1. Linear model on Average daily feed intake

\[y=X\beta + e\] Where:

\(y\) is a \(n x 1\) vector of Average feed intake (\(AFI_i\) kg) of an individual

\(X\beta\) are the fixed effects, \(location-trial \ (14)\)

\(e\) (\(n\) x \(1\)) is the random residuals vector.

and \(y_i-x_i\hat{\beta}=\hat{e_i}\) is the AFI of the animal \(i\)

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)
Estimated feed intake by location
Estimate Std. Error t value Pr(>|t|)
Loc_trialLoc_1_t_1 2.708833 0.0832692 32.53103 0
Loc_trialLoc_1_t_2 2.183523 0.0869718 25.10609 0
Loc_trialLoc_1_t_3 2.822786 0.0869718 32.45631 0
Loc_trialLoc_1_t_4 2.416517 0.1019835 23.69517 0
Loc_trialLoc_1_t_5 2.620605 0.1019835 25.69636 0
Loc_trialLoc_1_t_6 2.212822 0.0961510 23.01403 0
Loc_trialLoc_1_t_7 3.112632 0.0961510 32.37233 0
Loc_trialLoc_2_t_1 2.657176 0.0832692 31.91066 0
Loc_trialLoc_2_t_2 2.128764 0.0869718 24.47647 0
Loc_trialLoc_2_t_3 2.579544 0.0869718 29.65953 0
Loc_trialLoc_2_t_4 2.221091 0.1019835 21.77892 0
Loc_trialLoc_2_t_5 2.432934 0.1019835 23.85614 0
Loc_trialLoc_2_t_6 2.232878 0.0961510 23.22262 0
Loc_trialLoc_2_t_7 2.872766 0.1019835 28.16892 0 
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\\ less\ than\ or\ equal\ to\ 60 s\)

\(FollBLUP = animal\ ramdom\ effect\ as\ next\ individual\ that\ visited\ the\ feeder\\ less\ than\ or\ equal\ to\ 60 s\)

\(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\ estimated\ by\ random\ regression\)

\(Weight.gain.rr = group\ slope + trial.day.rr,\ estimated\ by\ random\ regression\)

\(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)
EarBLUP FollBLUP ADG ADG.locmean weight.gain.qr trial.day.rr Weight.gain.rr Average.daily.intake feed.intake.fitloc
EarBLUP 1.0000000 0.1109689 0.0760061 0.1387252 0.0342374 0.0956561 0.0580141 0.0637901 0.0929729
FollBLUP 0.1109689 1.0000000 -0.0044945 -0.0082033 0.0077479 -0.0017983 -0.0010907 0.0036300 0.0052907
ADG 0.0760061 -0.0044945 1.0000000 0.5478900 0.9250629 0.4983372 0.9567826 0.7109197 0.3885109
ADG.locmean 0.1387252 -0.0082033 0.5478900 1.0000000 0.5723908 0.9095570 0.5516344 0.4865273 0.7091039
weight.gain.qr 0.0342374 0.0077479 0.9250629 0.5723908 1.0000000 0.6489043 0.9870484 0.7318058 0.4860840
trial.day.rr 0.0956561 -0.0017983 0.4983372 0.9095570 0.6489043 1.0000000 0.6064868 0.5149118 0.7504736
Weight.gain.rr 0.0580141 -0.0010907 0.9567826 0.5516344 0.9870484 0.6064868 1.0000000 0.7302546 0.4551523
Average.daily.intake 0.0637901 0.0036300 0.7109197 0.4865273 0.7318058 0.5149118 0.7302546 1.0000000 0.6861158
feed.intake.fitloc 0.0929729 0.0052907 0.3885109 0.7091039 0.4860840 0.7504736 0.4551523 0.6861158 1.0000000

7. Visit length by animal

7.1. linear model of Visit length mean of each individual fitted by group mean

The sum of visit length total at the feeder for each individual was calculated as:

\[Vl_i= \sum_{j=1}^{n_i}X_j\]

Where:

\(Vl_i\) is the sum of the visit length total time at the feeder from the \(ith\) animal

\(X_j\) is the \(jth\) observation of visit length time at the feeder from the \(ith\) animal

\(n_i\) total number of observations from the \(ith\) animal

Visit length total time as deviation of estimated mean of location-trial group

The linear model fitted was: \[y=X\beta + e\] Where:

\(y\) is a \(n x 1\) vector of the sum of the visit length total \(Vl_i\) at the feeder from the \(ith\) animal

\(X\beta\) are the fixed effects, \(location-trial \ (12)\)

\(e\) (\(n\) x \(1\)) is the random residuals vector.

and \(y_i-x_i\hat{\beta}=\hat{e_i}\) is the visit length total as deviation of estimated mean of location-trial group of the animal \(i\)

visitID.mean<-trials.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)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1588    3159    5151    4968    6398    9643
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)
Mean estimated by location-trial group
Estimate Std. Error t value Pr(>|t|)
Loc_trialLoc_1_t_1 2056.836 292.2853 7.037085 0
Loc_trialLoc_1_t_2 5481.055 305.2820 17.954071 0
Loc_trialLoc_1_t_3 6575.724 305.2820 21.539837 0
Loc_trialLoc_1_t_4 4858.150 357.9749 13.571204 0
Loc_trialLoc_1_t_5 3241.925 357.9749 9.056292 0
Loc_trialLoc_1_t_6 6372.876 337.5019 18.882487 0
Loc_trialLoc_1_t_7 5843.033 337.5019 17.312592 0
Loc_trialLoc_2_t_1 2169.504 292.2853 7.422558 0
Loc_trialLoc_2_t_2 5314.229 305.2820 17.407607 0
Loc_trialLoc_2_t_3 7418.902 305.2820 24.301799 0
Loc_trialLoc_2_t_4 5624.471 357.9749 15.711916 0
Loc_trialLoc_2_t_5 2991.698 357.9749 8.357285 0
Loc_trialLoc_2_t_6 7056.528 337.5019 20.908110 0
Loc_trialLoc_2_t_7 5016.675 357.9749 14.014043 0 
visitID.mean$Loc_trial<-as.character(visitID.mean$Loc_trial)

#fitted by mean location
visitl.fitloc<-as.data.frame(vlm1$residuals)%>%
  mutate(., visitl.fitloc=vlm1$residuals,
         Eartag_trial=rownames(.))%>%select(Eartag_trial,visitl.fitloc)
ADGs.FI.BLUPs<-left_join(ADGs.FI.BLUPs, visitID.mean, by=c("Eartag_trial","Loc_trial"))%>%
  left_join(visitl.fitloc,by="Eartag_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<-trials.data%>%group_by(Loc_trial,Eartag_trial, trial.day)%>%
  summarize(visit.day=sum(visit.length))%>%ungroup()
vlm.day.lm1<-lmer(visit.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| = 2.20598 (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)
Estimates of fixed effects
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 100.3018761 6.6650317 251.8958 15.0489720 0.0000000
Loc_trialLoc_1_t_2 -7.6722271 9.2105712 210.3701 -0.8329806 0.4058004
Loc_trialLoc_1_t_3 10.0170791 9.2050543 209.8403 1.0882151 0.2777484
Loc_trialLoc_1_t_4 5.5997594 10.0514716 208.7188 0.5571084 0.5780503
Loc_trialLoc_1_t_5 8.3037668 10.3073719 230.6763 0.8056144 0.4212949
Loc_trialLoc_1_t_6 5.6554247 9.6588666 204.3212 0.5855164 0.5588467
Loc_trialLoc_1_t_7 54.4265445 9.7792982 214.6562 5.5654857 0.0000001
Loc_trialLoc_2_t_1 -2.8587496 9.4285765 252.1915 -0.3032006 0.7619870
Loc_trialLoc_2_t_2 -11.8623821 9.2102546 210.3414 -1.2879538 0.1991772
Loc_trialLoc_2_t_3 17.2303586 9.1918855 208.6758 1.8745184 0.0622562
Loc_trialLoc_2_t_4 11.7874008 10.0481893 208.4475 1.1730871 0.2420997
Loc_trialLoc_2_t_5 11.8787315 10.3073719 230.6763 1.1524501 0.2503291
Loc_trialLoc_2_t_6 25.0482761 9.6579622 204.2451 2.5935364 0.0101870
Loc_trialLoc_2_t_7 4.5536606 10.1002609 212.7824 0.4508458 0.6525597
trial.day -1.0337949 0.1851867 532.6249 -5.5824455 0.0000000
Loc_trialLoc_1_t_2:trial.day 0.8932764 0.2238279 264.3021 3.9909075 0.0000853
Loc_trialLoc_1_t_3:trial.day 0.8329100 0.2236384 263.1154 3.7243607 0.0002396
Loc_trialLoc_1_t_4:trial.day 0.6067473 0.2397692 243.6376 2.5305469 0.0120193
Loc_trialLoc_1_t_5:trial.day 0.6018994 0.2634306 353.1898 2.2848504 0.0229140
Loc_trialLoc_1_t_6:trial.day 0.6649842 0.2302288 237.8110 2.8883626 0.0042297
Loc_trialLoc_1_t_7:trial.day 0.0602560 0.2375870 269.3712 0.2536166 0.7999852
Loc_trialLoc_2_t_1:trial.day 0.3607581 0.2622868 535.4997 1.3754340 0.1695720
Loc_trialLoc_2_t_2:trial.day 0.9290088 0.2238117 264.2269 4.1508507 0.0000447
Loc_trialLoc_2_t_3:trial.day 0.7954158 0.2229391 260.1725 3.5678619 0.0004284
Loc_trialLoc_2_t_4:trial.day 0.8223142 0.2395548 242.7746 3.4326768 0.0007025
Loc_trialLoc_2_t_5:trial.day 0.1369993 0.2634306 353.1898 0.5200583 0.6033487
Loc_trialLoc_2_t_6:trial.day 0.4300173 0.2302195 237.7731 1.8678577 0.0630121
Loc_trialLoc_2_t_7:trial.day 1.0945119 0.2432637 258.0200 4.4992812 0.0000103 
# Variances
print(VarCorr(vlm.day.lm1), comp="Variance")
##  Groups       Name        Variance  Corr  
##  Eartag_trial (Intercept) 411.22665       
##               trial.day     0.15582 -0.566
##  Residual                 373.97377
# fixed effects
fxeff<-as.matrix(summary(vlm.day.lm1)$coef[,1])
fxeff<-tibble(Loc_trial=rownames(fxeff)[1:14],
              intercept.vl=fxeff[1:14,],slope.vl=fxeff[15:28,])%>%
        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%>%kable()%>%kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F, position = "left", font_size = 12)
Loc_trial intercept.vl slope.vl
Loc_1_t_1 100.30188 -1.0337949
Loc_1_t_2 92.62965 -0.1405185
Loc_1_t_3 110.31896 -0.2008848
Loc_1_t_4 105.90164 -0.4270476
Loc_1_t_5 108.60564 -0.4318955
Loc_1_t_6 105.95730 -0.3688107
Loc_1_t_7 154.72842 -0.9735389
Loc_2_t_1 97.44313 -0.6730367
Loc_2_t_2 88.43949 -0.1047861
Loc_2_t_3 117.53223 -0.2383791
Loc_2_t_4 112.08928 -0.2114807
Loc_2_t_5 112.18061 -0.8967956
Loc_2_t_6 125.35015 -0.6037776
Loc_2_t_7 104.85554 0.0607171 
# 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))
Estimates correlation between BLUP’s,visit length, weight gain and feed intake
Eartag Blup
Follower Blup
Correlation p-val Correlation p-val
trial.day.rr 0.0956561 0.2697576 -0.0017983 0.9834848
feed.intake.fitloc 0.0929729 0.2834835 0.0052907 0.9514380
visitl.fitloc 0.4408095 0.0000001 -0.2227322 0.0094164
Visit.length -0.0215532 0.8040371 -0.1374602 0.1118692
vl.day -0.0299928 0.7298523 -0.1912855 0.0262541 
#output.vl%>%select("EarBLUP","FollBLUP","trial.day.rr","feed.intake.fitloc",
 #                  "visitl.fitloc", "Visit.length","vl.day")%>%cor()