Clumpsy code for calculating mean and standard error

Does this give you what you are looking for? It seems to reproduce the output. The note there is that I think your result actually has the columns labeled wrong. One of your steps did a character sort on the columns, so they are not in the order you expected when your renamed later.

What we do is first use pivot_longer() so that we have tidy data. From there, we can group and summarize to calculate the mean and standard errors. Then we pivot_wider() to move it back to the wide format of your result.

library(dplyr)
library(tidyr)
library(plotrix)

col_order <- paste0(rep(paste0(c("mean", "se"), "_day_"), length(unique(df1$days_incubated))),
                    rep(unique(df1$days_incubated), each = 2))

df1 %>% 
  pivot_longer(-days_incubated, names_to = "mass") %>% 
  group_by(days_incubated, mass) %>% 
  summarize(mean = mean(value),
            se = std.error(value), .groups = "drop") %>%
  pivot_wider(names_from = days_incubated, values_from = c("mean", "se"),
              names_glue = "{.value}_day_{days_incubated}") %>% 
  relocate(mass, all_of(col_order))

# A tibble: 5 x 19
  mass  mean_day_0 se_day_0 mean_day_4 se_day_4 mean_day_10 se_day_10 mean_day_17 se_day_17 mean_day_24 se_day_24 mean_day_66 se_day_66 mean_day_81
  <chr>      <dbl>    <dbl>      <dbl>    <dbl>       <dbl>     <dbl>       <dbl>     <dbl>       <dbl>     <dbl>       <dbl>     <dbl>       <dbl>
1 i.x3~     24.2    13.9        45.1    17.5          3.29    0.932        26.6     19.5        113.       68.6        10.1       4.67       13.0  
2 i.x3~      0.036   0.0287      0.048   0.0116       0.014   0.00748       0.06     0.0429       0.198     0.121       0.016     0.006       0.01 
3 i.x3~      0.178   0.0867      0       0            0       0             0        0            0.01      0.01        0         0           0    
4 i.x3~      0.02    0.0155      0       0            0       0             0        0            0.008     0.008       0.002     0.002       0.006
5 i.x4~      0.182   0.0978      0.294   0.240        0.174   0.0941        0.158    0.119        0.378     0.210       0.006     0.004       0.018
# ... with 5 more variables: se_day_81 <dbl>, mean_day_94 <dbl>, se_day_94 <dbl>, mean_day_116 <dbl>, se_day_116 <dbl>

Line-by-line

First, I have added in a short step to help column sorting.

• col_order is a character vector we create of the correct final order we are looking for. Basically we are just concatenating strings to match the column names. You can skip this and simply type each name individually much like you did in your select() statement, but this saves typing.

Now onto the pipeline.

• pivot_longer() will take the data from the columns and stack them. The column names will now be in a variable called mass. By default, the values will be in a new variable called value.
• group_by() groups on the specified variables, which are the row and column combinations you want in the final table.
• summarize() collapses the multiple rows for each mass and days_incubated combination into one new row, with two columns: the mean and se. The .groups = "drop" argument says to drop the groups (i.e., the tibble is no longer grouped, check the help file for more details).
• pivot_wider() takes the long formatted summary tibble and pivots the data out into columns, much how you had your format. It says to take the column names from days_incubated, and to take the values from the mean and the se columns. The names_glue argument specifies how to name the new columns. It uses syntax from the glue packages, there the { brackets basically mean "substitute a value here". So it is the '[name of value column]_day_[days_incubated]' as each name.
• relocate() simply reorders the columns. This is relatively new to dplyr, and is the preferred way now to rearrange columns instead of select(). It says to first take mass, and then all_of() the columns in the col_order vector that we created earlier.

I'd like to offer a solution mainly in base R.

library(dplyr)

# Define a function that computes the standard error of a vector.
f.SE <- function(x) sd(x)/sqrt(length(x))

# Define a function that compute the standard error of a column (MARGIN = 2 for column). 
f.colSE <- function(mydf) apply(mydf, MARGIN = 2, FUN = f.SE)

# Define a function that combines the means and the standard errors of a column
f.colMeanSE <- function(mycol) cbind(means = colMeans(mycol), se = f.colSE(mycol))

# Apply the function to all numeric columns in your data frame by using the days_incubated subgroup as the indices. The result is a list. 
my_mean_se1 <- by(data = df1[,-1], INDICES = df1$days_incubated, FUN= f.colMeanSE); my_mean_se1

#df1$days_incubated: 0
#means          se
#i.x33.031 24.182 13.91245823
#i.x35.034  0.036  0.02874022
#i.x36.017  0.178  0.08668333
#i.x39.959  0.020  0.01549193
#i.x40.023  0.182  0.09779571
#------------------------------------------------------------ 
#  df1$days_incubated: 10
#means          se
#i.x33.031 3.294 0.932119091
#i.x35.034 0.014 0.007483315
#i.x36.017 0.000 0.000000000
#i.x39.959 0.000 0.000000000
#i.x40.023 0.174 0.094106323
#------------------------------------------------------------

# Bind the columns. 
my_mean_se2 <- do.call(cbind, my_mean_se1); my_mean_se2

#           means          se means          se means          se  means          se   means         se  means          se  means       se  means
#i.x33.031 24.182 13.91245823 3.294 0.932119091 0.036 0.036000000 26.628 19.46476571 113.340 68.5767486 45.140 17.49416646 10.144 4.665935 13.036
#i.x35.034  0.036  0.02874022 0.014 0.007483315 0.008 0.004898979  0.060  0.04289522   0.198  0.1212601  0.048  0.01157584  0.016 0.006000  0.010
#i.x36.017  0.178  0.08668333 0.000 0.000000000 0.028 0.023323808  0.000  0.00000000   0.010  0.0100000  0.000  0.00000000  0.000 0.000000  0.000
#i.x39.959  0.020  0.01549193 0.000 0.000000000 0.006 0.004000000  0.000  0.00000000   0.008  0.0080000  0.000  0.00000000  0.002 0.002000  0.006
#i.x40.023  0.182  0.09779571 0.174 0.094106323 0.018 0.015620499  0.158  0.11939012   0.378  0.2103188  0.294  0.23972067  0.006 0.004000  0.018
 #                  se means         se
#i.x33.031 6.002108463 0.094 0.07304793
#i.x35.034 0.007745967 0.002 0.00200000
#i.x36.017 0.000000000 0.136 0.07194442
#i.x39.959 0.004000000 0.024 0.01122497
#i.x40.023 0.011135529 0.016 0.00600000

# Bind the names of the columns to the names of days_incubated
names_days <- unique(df1$days_incubated)
names_meanse <- rbind(paste0("mean day ", names_days), paste0("se day ", names_days))
colnames(my_mean_se2) <- names_meanse 

# Round the numbers to 2 decimal places and convert the numeric table to a data frame.
mean_se <- as.data.frame(round(my_mean_se2, digits = 2)); mean_se

#         mean day 0 se day 0 mean day 4 se day 4 mean day 10 se day 10 mean day 17
#i.x33.031      24.18    13.91       3.29     0.93        0.04      0.04       26.63
#i.x35.034       0.04     0.03       0.01     0.01        0.01      0.00        0.06
#i.x36.017       0.18     0.09       0.00     0.00        0.03      0.02        0.00
#i.x39.959       0.02     0.02       0.00     0.00        0.01      0.00        0.00
#i.x40.023       0.18     0.10       0.17     0.09        0.02      0.02        0.16
#          se day 17 mean day 24 se day 24 mean day 66 se day 66 mean day 81 se day 81
#i.x33.031     19.46      113.34     68.58       45.14     17.49       10.14      4.67
#i.x35.034      0.04        0.20      0.12        0.05      0.01        0.02      0.01
#i.x36.017      0.00        0.01      0.01        0.00      0.00        0.00      0.00
#i.x39.959      0.00        0.01      0.01        0.00      0.00        0.00      0.00
#i.x40.023      0.12        0.38      0.21        0.29      0.24        0.01      0.00
#          mean day 94 se day 94 mean day 116 se day 116
#i.x33.031       13.04      6.00         0.09       0.07
#i.x35.034        0.01      0.01         0.00       0.00
#i.x36.017        0.00      0.00         0.14       0.07
#i.x39.959        0.01      0.00         0.02       0.01
#i.x40.023        0.02      0.01         0.02       0.01

# Convert the data frame to a tibble with rowname set as a new column
as_tibble(mean_se, rownames = "rowname")

# A tibble: 5 x 19
#  rowname `mean day 0` `se day 0` `mean day 4` `se day 4` `mean day 10` `se day 10` `mean day 17` `se day 17` `mean day 24` `se day 24`
#  <chr>          <dbl>      <dbl>        <dbl>      <dbl>         <dbl>       <dbl>         <dbl>       <dbl>         <dbl>       <dbl>
#1 i.x33.~        24.2       13.9          3.29       0.93          0.04        0.04         26.6        19.5         113.         68.6 
#2 i.x35.~         0.04       0.03         0.01       0.01          0.01        0             0.06        0.04          0.2         0.12
#3 i.x36.~         0.18       0.09         0          0             0.03        0.02          0           0             0.01        0.01
#4 i.x39.~         0.02       0.02         0          0             0.01        0             0           0             0.01        0.01
#5 i.x40.~         0.18       0.1          0.17       0.09          0.02        0.02          0.16        0.12          0.38        0.21
# ... with 8 more variables: `mean day 66` <dbl>, `se day 66` <dbl>, `mean day 81` <dbl>, `se day 81` <dbl>, `mean day 94` <dbl>, `se day
#   94` <dbl>, `mean day 116` <dbl>, `se day 116` <dbl>

I think what you started off with was in right direction. You can next get data in long format, change the column names the way you want it and get the data back in wide format.

library(dplyr)
library(tidyr)

df1 %>%
  group_by(days_incubated = as.numeric(days_incubated)) %>%
  summarise(across(i.x33.031:i.x40.023,  list(mean = mean, se = plotrix::std.error))) %>%
  pivot_longer(cols = -days_incubated, 
               names_to = c('mass', 'col'), 
               names_sep = '_') %>%
  #If you need column exactly as shown
  mutate(col = paste(col, 'day')) %>%
  pivot_wider(names_from = c(col, days_incubated), values_from = value, names_sep = ' ')


#  mass  `mean day 0` `se day 0` `mean day 4` `se day 4` `mean day 10` `se day 10`
#  <chr>        <dbl>      <dbl>        <dbl>      <dbl>         <dbl>       <dbl>
#1 i.x3…       24.2      13.9          45.1      17.5            3.29      0.932  
#2 i.x3…        0.036     0.0287        0.048     0.0116         0.014     0.00748
#3 i.x3…        0.178     0.0867        0         0              0         0      
#4 i.x3…        0.02      0.0155        0         0              0         0      
#5 i.x4…        0.182     0.0978        0.294     0.240          0.174     0.0941 
# … with 12 more variables: `mean day 17` <dbl>, `se day 17` <dbl>, `mean day
#   24` <dbl>, `se day 24` <dbl>, `mean day 66` <dbl>, `se day 66` <dbl>, `mean day
#   81` <dbl>, `se day 81` <dbl>, `mean day 94` <dbl>, `se day 94` <dbl>, `mean day
#   116` <dbl>, `se day 116` <dbl>