Create lm's from split data

Question

Essentially what I want to do is automate this for all my data:

plants_A<-plants_sorted[plants_sorted[, 'treatment']== 'A', ]
plants_A1<-plants_A[plants_A[, 'replicate']== '1', ]
lm(weight~time, data = plants_A1)

From 'plants' I want to make lm's for all treatment and replicate combinations.

I have also managed to split the data using:

plants_treat_repl <- split(plants, paste(plants$treatment, plants$replicate))

But I can't seem to make lm's from this split data.

# sample data
structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 1L, 2L ), .Label = c("A", 
"B", "C", "D"), class = "factor"), replicate = c(1, 2, 3, 4, 5, 1), time = 
structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Label = c("6", "8", "10", "12", "14", 
"16"), class = "factor"), weight = c(2, 0, 0, 0, 0.5, 2.6), trtrep = 
structure(c(1L, 5L, 9L, 13L, 17L, 2L), .Label = c("A.1", "B.1", "C.1", "D.1", 
"A.2", "B.2", "C.2", "D.2", "A.3", "B.3", "C.3", "D.3", "A.4", "B.4", "C.4", 
"D.4", "A.5", "B.5", "C.5", "D.5"), class = "factor")), row.names = c(NA, 6L), 
class = "data.frame")

I have 4 treatments: A-D and 5 replicates for each treatment. I want to make models for each combination of treatment/replicate: A1, A2, ... D5. — Martijn Ekhart, Mar 12 '21 at 14:02
Try `list_of_lms = lapply(plants_treat_repl, lm, formula = weight ~ time)` — Gregor Thomas, Mar 12 '21 at 14:06
@GregorThomas Thank you for the response! That also seems to work and give similar results to Ben's code but it results in the same problem. I still only get a list of values and not linear models out of it. (could also be my lack of statistical background that I am missing something) — Martijn Ekhart, Mar 12 '21 at 14:16
No, that returns a list of linear models. Check with `lapply(list_of_lms, class)`. If you want the coefficients, `lapply(list_of_lms, coef)`. — Gregor Thomas, Mar 12 '21 at 14:27
@GregorThomas I still get the same thing back as coefficients as I get from other methods. It still just shows the weights at the different times and not slopes/intercepts. — Martijn Ekhart, Mar 12 '21 at 14:35
I've added a fully reproducible answer based on the built-in `mtcars` data set. If you need more help figuring out how to get things working on your data, please share a reproducible sample of data -- using `dput()` to make it copy/pasteable is the nicest way to share sample data. — Gregor Thomas, Mar 12 '21 at 14:39
I edited your sample data into your question. In the sample data you show, a treatment:replicate interaction is a single row. So you're trying to fit a separate model to each individual row of data - this will not work. You need multiple observations (rows) to fit a model. — Gregor Thomas, Mar 12 '21 at 15:12
Sorry, I am very confused now. When I do it manually with 'plants_A1' I get a table with 6 rows: time/weight for A1 at 6 times. Then with (manually) with --- lm_A1<-lm(weight~time, data = plants_A1) --- I still don't get a lm out of it, even with multiple rows. — Martijn Ekhart, Mar 12 '21 at 15:17
data plants_A1: --- structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 1L, 1L ), .Label = c("A", "B", "C", "D"), class = "factor"), replicate = c(1, 1, 1, 1, 1, 1), time = structure(1:6, .Label = c("6", "8", "10", "12", "14", "16"), class = "factor"), weight = c(2, 7.9, 17.8, 31.5, 44.5, 59.9)), row.names = c(NA, -6L), class = c("tbl_df", "tbl", "data.frame")) --- — Martijn Ekhart, Mar 12 '21 at 15:19
I am getting even more confused, when knitting my document in R markdown, the PDF does show my correct coefficients but my console inside Rstudio does not. — Martijn Ekhart, Mar 12 '21 at 15:36
I think the issue may be that your `time` column is `factor` (categorical) when it should be `numeric`. In your original data, convert it to numeric with `plants_sorted$time = as.numeric(as.character(plants_sorted$time))` and then try again. It's possible you accidentally converted the column at some point in your RStudio session, but did not put the conversion code in your knitr document. — Gregor Thomas, Mar 12 '21 at 15:37
Ah, that should be it. In some other exercises I had to consider it as a factor but should have changed it back. Why it is different in the pdf, I still wouldn't know. Thank you for the help! — Martijn Ekhart, Mar 12 '21 at 16:11

score 1 · Answer 1 · answered Mar 12 '21 at 14:19

1

Try lapply function:

lms <- lapply(split(plants, paste(plants$treatment,plants$replicate)),
function(x)lm(weight~time, data = x))

lms is a list of all posibles lm's.

answered Mar 12 '21 at 14:19

Marcos Pérez

1,260
2
7

Thank you for the answer! This seems to work very well except I still don't fully get what I was after. When checking it with 'class(lms)' it also says 'lm' but I can't seem to get a value for slope out of it. I usually would do this with 'coef(lms)' but this just returns the values of the weight at the specific times and not a slope/intercept etc. – Martijn Ekhart Mar 12 '21 at 14:25
You can do it with `lapply(lms,coef)`. If you want to predict you can use `lapply(lms,function(lm)predict(lm,new_data))`. – Marcos Pérez Mar 12 '21 at 14:55

score 1 · Accepted Answer · answered Mar 12 '21 at 14:38

1

Let's transition this to a reproducible example using built-in data. (Feel free to share reproducible sample data in the question - dput(your_data[1:10, ]) is a great way to share data reproducibly, pick an appropriate subset.)

## split the data by group
mt_split = split(mtcars, mtcars$cyl)

## fit models to each group
mods = lapply(mt_split, lm, formula = mpg ~ wt)

## extract the coefficients from each model
lapply(mods, coef)
# $`4`
# (Intercept)          wt 
#   39.571196   -5.647025 
# 
# $`6`
# (Intercept)          wt 
#   28.408845   -2.780106 
# 
# $`8`
# (Intercept)          wt 
#   23.868029   -2.192438

answered Mar 12 '21 at 14:38

Gregor Thomas

136,190
20
167
294

From the car package I get correct results but when substituting for my data (also when just selecting for treatment, not treatment and replicate) I still only get my values back and no intercepts/slopes. There is probably something wrong with my data then. If you know what could possibly be wrong with it then feel free to say so but otherwise still thanks so much for the help! – Martijn Ekhart Mar 12 '21 at 14:45
1

I'd be happy to help diagnose your data issues - you just need to share some of your data reproducibly as I describe in my answer and in the comments on your question. If you need more help making a reproducible example, [this FAQ goes into great detail](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). – Gregor Thomas Mar 12 '21 at 14:51
I have tried to do it but it gives errors with that the code is not formatted correctly and since my deadline is in about an hour I don't have time to get into that unfortunately. I will get back to you after that if you would still want to help because I am getting curious. – Martijn Ekhart Mar 12 '21 at 15:04
structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 1L, 2L ), .Label = c("A", "B", "C", "D"), class = "factor"), replicate = c(1, 2, 3, 4, 5, 1), time = structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Label = c("6", "8", "10", "12", "14", "16"), class = "factor"), weight = c(2, 0, 0, 0, 0.5, 2.6), trtrep = structure(c(1L, 5L, 9L, 13L, 17L, 2L), .Label = c("A.1", "B.1", "C.1", "D.1", "A.2", "B.2", "C.2", "D.2", "A.3", "B.3", "C.3", "D.3", "A.4", "B.4", "C.4", "D.4", "A.5", "B.5", "C.5", "D.5"), class = "factor")), row.names = c(NA, 6L), class = "data.frame") – Martijn Ekhart Mar 12 '21 at 15:05
That is the shortest I can get here, maybe that works? When putting it in the original post I get the error message. – Martijn Ekhart Mar 12 '21 at 15:05

score 0 · Answer 3 · answered Mar 12 '21 at 14:05

0

You haven't told us exactly what doesn't work, but:

plants <- transform(plants, trtrep = interaction(treatment, replicate)
lme4::lmList(weight~time | trtrep, data=plants)

should work. You might get away with defining the replicate on the fly:

lme4::lmList(weight~time | interaction(treatment, replicate), data=plants)

answered Mar 12 '21 at 14:05

Ben Bolker

211,554
25
370
453

Thanks a lot for the swift response! I now do get a df but this only contains the intercepts and weights at the different time values. How would I get slopes/equations out of this? – Martijn Ekhart Mar 12 '21 at 14:13

Create lm's from split data

3 Answers3