Linear regression on subsets of data corresponding to two factor combination

Question

[botinky]

In R I am trying to execute linear regressions of Predator.mass/Prey.mass within the subsets of Lifestage and Feeding Interaction. The data I have currently look like this:

|      | Predator.lifestage | Type.of.feeding.interaction | Predator.mass | Prey.mass |
|------|--------------------|-----------------------------|---------------|-----------|
| 1    | adult              | predacious/piscivorous      | 1540          | 14.3      |
| 2    | adult              | predacious/piscivorous      | 1600          | 6.02      |
| 3    | adult              | predacious/piscivorous      | 1840          | 11.9      |
| 4    | adult              | predacious/piscivorous      | 87.6          | 8.12      |
| 5    | adult              | predacious/piscivorous      | 63.9          | 6.56      |
| 6    | adult              | predacious/piscivorous      | 79.2          | 5.41      |
| 7    | adult              | predacious/piscivorous      | 71.2          | 4.45      |
| 8    | adult              | predacious/piscivorous      | 92.1          | 5.98      |

The list of lifestages goes down to include a total of six factors (adult, juvenline etc.) and feeding interactions go down to include a total of five factors (predacious/piscivorous, insectivorous etc.).

So what I want to achieve is a linear regression of Predator.mass/Prey.mass for each combination of subsets (adult, predacious/piscivorous; juvenile, predacious/piscivorous; lavar, predacious/piscivorous... and so on through the lifestages and feeding interactions).

What I'm trying to do is use a for loop and this is what I've got so far...

for (Type.of.feeding.interaction in My_Data) {
                for (Predator.lifestage in My_Data) {
                LinMod <- lm(Predator.mass ~ Prey.mass)
        }
      }

When I run this says object 'Predator.mass' not found. I've only been using coding a couple of weeks so I'm still fumbling around in the dark a bit. I'd really appreciate any help with this! Thank you :)

[cpad0112]

try :

with (My_Data, for (Type.of.feeding.interaction in My_Data) {
                for (Predator.lifestage in My_Data) {
                LinMod <- lm(Predator.mass ~ Prey.mass)
        }
      }
)

is this what you are looking for:

> df
  Predator.lifestage Type.of.feeding.interaction Predator.mass Prey.mass
1              adult      predacious/piscivorous        1540.0     14.30
2              adult      predacious/piscivorous        1600.0      6.02
3              adult      predacious/piscivorous        1840.0     11.90
4              adult      predacious/piscivorous          87.6      8.12
5              adult      predacious/piscivorous          63.9      6.56
6              adult      predacious/piscivorous          79.2      5.41
7              adult      predacious/piscivorous          71.2      4.45
8              adult      predacious/piscivorous          92.1      5.98

> library(dplyr)
> library(broom)
> df %>%
       group_by(Type.of.feeding.interaction, Predator.lifestage) %>%
       group_map(~ tidy(lm(Predator.mass ~ Prey.mass, .)))
# A tibble: 2 x 7
# Groups:   Type.of.feeding.interaction, Predator.lifestage [1]
  Type.of.feeding.interaction Predator.lifestage term        estimate std.error statistic p.value
  <chr>                       <chr>              <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 predacious/piscivorous      adult              (Intercept)    -618.     593.      -1.04  0.338 
2 predacious/piscivorous      adult              Prey.mass       164.      69.9      2.35  0.0568

[botinky]

Hi there, cheers for your reply! I've applied it to my data set and it works, just as you said. Thanks again and have a great day :)

Answer 1

[Sam]

Assuming your file is called data.txt

library(data.table)
dat <- fread("data.txt")
results <- dat[,.(Pvalue=summary(lm(Predator.mass~Prey.mass))$coefficient[2,4], 
          beta=summary(lm(Predator.mass~Prey.mass))$coefficient[2,3], 
          r2=summary(lm(Predator.mass~Prey.mass))$r.squared),
          by=c("Predator.lifestage","Type.of.feeding.interaction")]

Should give you the expected results.

[botinky]

Hi Sam, thank you so much for taking the time to respond to my question! I'm affraid when I tried to run the code I got this error message: Error in summary(lm(Predator.mass ~ Prey.mass))$coefficient[2, 4] : subscript out of bounds