Well, as is often said in R tutorials, "if you're using a for loop, you're doing it wrong." In addition to leaving the code a little slower, you are not using one of the most important factors of language, which is its functional orientation, extremely useful in problems like yours: apply complex functions structurally into data.

Robert's answer is correct. I only offered a more modern version using the dplyr + tidyr , broom and purrr packages to:

The structure of my code should help you to estimate your models without problems.

library(tidyverse)

## gerando dados com a mesma estrutura apresentada

set.seed(134)

df <- 
  tibble(Animal = 1:190, taxa = runif(n = 190, min = 3, max = 23)) %>%  
  crossing(dia = 0:3) %>%
  mutate(consumo = 245 + dia*taxa)


df


# A tibble: 760 x 4
   Animal    taxa   dia  consumo
   <int>     <dbl> <int>    <dbl>
1      1  7.007272     0 245.0000
2      1  7.007272     1 252.0073
3      1  7.007272     2 259.0145
4      1  7.007272     3 266.0218
5      2 15.551850     0 245.0000
6      2 15.551850     1 260.5518
7      2 15.551850     2 276.1037
8      2 15.551850     3 291.6555
9      3 20.629976     0 245.0000
10     3 20.629976     1 265.6300
# ... with 750 more rows

We now estimate a simple linear model with lm for the consumo ~ dia function and extract the coefficients of each model with the help of the broom::tidy function, which extracts the coefficients and statistics in a data.frame "tidy" (clean):

## gerando um modelo linear no dia para cada animal e extraindo os coeficientes  --------------

df %>% 
  group_by(Animal) %>% 
  nest(.key = dados_animal) %>% 
  mutate(model = map(dados_animal, ~lm(consumo ~ dia, data = .))) %>% 
  mutate(coef = map(model, broom::tidy)) %>% 
  unnest(coef)

Resulting in:

# A tibble: 380 x 6
   Animal        term   estimate    std.error    statistic      p.value
   <int>       <chr>      <dbl>        <dbl>        <dbl>        <dbl>
1      1 (Intercept) 245.000000 1.174950e-15 2.085196e+17 2.299886e-35
2      1         dia   7.007272 6.280370e-16 1.115742e+16 8.032903e-33
3      2 (Intercept) 245.000000 0.000000e+00          Inf 0.000000e+00
4      2         dia  15.551850 0.000000e+00          Inf 0.000000e+00
5      3 (Intercept) 245.000000 4.699798e-15 5.212990e+16 3.679818e-34
6      3         dia  20.629976 2.512148e-15 8.212086e+15 1.482836e-32
7      4 (Intercept) 245.000000 0.000000e+00          Inf 0.000000e+00
8      4         dia   8.025121 0.000000e+00          Inf 0.000000e+00
9      5 (Intercept) 245.000000 0.000000e+00          Inf 0.000000e+00
10     5         dia  15.198999 0.000000e+00          Inf 0.000000e+00
# ... with 370 more rows

To learn more about the "tidyverse" way of programming, I recommend taking a look at articles in the pacency site and in the book R for Data Sceince , especially in the part about working with many models , because that's basically how I presented this answer.

I invented data and a function, but can use the same logic:

set.seed(123)
table1 = data.frame(Animal=rep(1:6,each=10),Dia=rep(1:4,15), Consumo=round(runif(60,240,500),0))

runDLM = function(beta,dat) { 
  x=dat$Dia 
  y=dat$Consumo
  coef(lm(y~I(x^beta)))
  }

sdat=split(table1,table1$Animal)
betat=2
sapply(sdat,runDLM,beta=betat)


#                         1          2          3          4          5          6
#    (Intercept) 355.804959 319.156589 376.897521 398.223256 313.494215 369.547287
#    I(x^beta)     5.322314   6.699225   3.538843  -2.167442   2.523967  -1.993798