0

I tried to create a graph in R to describe a quadratic term with a random effect.

I used the function lme() with a random effect which included random slope and intercept:

    diversity = c(0.69, 1.54, 0.84, 1.48, 1.71, 1.80, 2.09, 1.63, 2.40, 2.20, 
        2.56, 2.30, 2.67, 1.98, 1.65, 2.33, 2.17, 1.98, 1.96, 1.33, 2.55, 2.49, 2.39, 2.47, 2.42, 2.44, 2.35, 2.33, 2.01, 2.39)
    Plot_age = c(7, 7, 9, 12, 17, 19, 22, 32, 31, 35, 35, 36, 36, 36, 37, 37,
 37, 38, 38, 38, 110, 111, 112, 113, 115, 116, 117, 118, 120, 121)
    subject = c("A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C")

this data frame gives you only an overview of my actual data. I hope it will work anyways.

age <- Plot_age
div <- diversity
age2 <- Plot_age^2
qm1r <- lme( fixed = div ~ age+age2, 
             random = ~ age |subject)

To get a graph I used the plot() function as usually. With the lm() it worked, now I also used the predict() function and lme() with the following terms:

plot(age, div, pch = 16)

X <- seq(0, 200, 0.1) #X has a range of 0 to 200 which include all of my x-axis data

NewData <- data.frame(age= Plot_age,
                      age2 = Plot_age^2,
                      subject = subject)

Y <- predict(qm1r, NewData)

points(Y ~ X, type ="l", lwd=3)

I got a graphical representation of my points() function. Unfortunately it didnt came out right (I wanted a quadratic parabola, no nonlinear line) -->see image below

xmax <- X[Y == max(Y)]

#The reason why I used this function was to get the xmax depending on the peak of the quadratic parabola graph (this worked when I only used the lm() function, lme() did not work) Is there any way to create a graphical representation of my lme() function with random effects? I would be so glad if someone could solve this problem :) I got a graph but no quadratic regression out of the points() function enter image description here

Gregor Thomas
  • 136,190
  • 20
  • 167
  • 294
Helene2022
  • 1
  • 1
  • *"also my newdata didnt work as I expected and I dont know why"* What did you expect, and how was it different from what you got? If you could share a but of sample data it would help a lot: something like `dput(Data1[1:20, ])` for the first 20 rows (which will include all class and structure information and be copy/pasteabled into R). – Gregor Thomas May 03 '22 at 21:03
  • 1
    As a general comment, your code does a lot of pulling vectors out of the data frame, things like `age <- Data1$Plot_age` creates `age` as its own vector. Why take it out of the data frame and give it a different name? It's just giving an opportunity for bugs – Gregor Thomas May 03 '22 at 21:04
  • @GregorThomas I am using the vector due to the fact that when I am starting a new dataframe for the predict() function there is always an error when I wrote Plot_age^2=X^2. Using the age2= Plot_age^2 there were no errors. – Helene2022 May 04 '22 at 11:07
  • The syntax of your sample data doesn't work. `subject = c(A, B, ...` needs quotes around all the letters, and it ends with a comma making in incomplete. If you use `dput()` and copy/paste the result into your question (it should start with `structure(list...`) the syntax will be right and the class information will be included. You can see the FAQ about [making reproducible examples in R](https://stackoverflow.com/a/5963610/903061) if you need more help with that. – Gregor Thomas May 04 '22 at 14:02
  • @GregorThomas oh yes, so sorry. I was just giving you an example. the subjects arent called like in this list. My problem is that I dont get a quadratic parabola and no xmax- values depending on x and y – Helene2022 May 04 '22 at 15:09
  • Yes, I understand the problem. It's just hard to demonstrate a solution without sample data to test on. Looks much better now--we can put it in a data frame and it's almost as easy as if you had use `dput()`. – Gregor Thomas May 04 '22 at 15:43

1 Answers1

0
## put your data in a data frame
df = data.frame(diversity, Plot_age, subject)    

## define the model using the data frame data
## note the `I()` which lets us square a variable in the formula
## so we can skip the step of creating a vector of squared values
library(nlme)
qm1r <- lme( fixed = diversity ~ Plot_age + I(Plot_age^2), 
             random = ~ Plot_age | subject,
             data = df)

## make the initial plot, again using data in the data frame
with(df, plot(Plot_age, diversity, pch = 16))

## create data for prediction, also in a data frame
## Use the same column name as original data so the model knows what to do!
## 201 points is plenty, setting `by = 0.1` to get 2001 points isn't necessary
pred_data = data.frame(
  Plot_age = seq(0, 200)
)

## add predictions to the pred_data data frame
## level = 0 is population-level
pred_data$diversity = predict(qm1r, newdata = pred_data, level = 0)

## add the predicted line to the plot
with(pred_data, lines(Plot_age, diversity))

enter image description here

If you want the x axis to go all the way to 200, set xlim = c(0, 200) in the initial plot() call.

Gregor Thomas
  • 136,190
  • 20
  • 167
  • 294
  • Wow, thank you so much! That helped me alot :) Last question: Do you know if I can get the xmax value of this predicted line depending on X and Y? – Helene2022 May 04 '22 at 20:38
  • Do you mean you want to know the x value corresponding to the maximum y value? The extreme value of a parabola ax^2 + bx + c is located at x = -b / 2a ([just googled it](https://www.google.com/search?client=firefox-b-1-d&q=how+to+find+max+of+parabola)). Using that notation `a` is the `Plot_age^2`, coefficient and `b` is the `Plot_age` coefficient, so have a look at `fixef(qm1r)` calculate the x-value of the max directly. You could then use `predict()` again to get the corresponding y value. – Gregor Thomas May 05 '22 at 00:43
  • If you consider your question resolved, please 'accept' the answer by clicking the checkmark next to it in the left margin. – Gregor Thomas May 05 '22 at 00:43