2
$\begingroup$

I would like to determine if there is a significant relationship between a measured continuous variable (predictor variable) and a response variable, which are both measured over time and therefore dependent on previous time points. Additionally, I conducted multiple experiments with different types of subjects, and I want to determine if the relationship exists across the subjects.

Simulated example: I want to measure the effect of soil reuse on plant growth (biomass). I conducted three experiments, each with a different plant species. The treatment group had soil plots with “reused” soil, which was reused consecutively in each round of the experiment to grow a new batch of plants. The control group had soil plots with “new” soil, which was new in each round of the experiment (so control soil plots were actually independent across rounds). Potassium concentrations in the treatment group soil plots were measured at the start of each round (control group had fixed potassium concentration). Plant biomass in each soil plot was measured at the end of each round, and this variable was used to calculate an effect size (specifically, the log-response ratio) that incorporates the control and treatment group measurements.

I want to determine if the predictor variable (initial potassium level in the reused soil) is significantly correlated with the response variable (biomass effect size), across multiple types of subjects (plant species).

What would be the appropriate strategy to determine if a significant relationship exists, considering the dependencies in the data?

Simulated example data in R:

set.seed(1000)
data <- data.frame(species = c( rep("species1", 5), rep("species2", 5), rep("species3", 5)), # species used in the experiment
                   round = rep(1:5,3), # rounds (equivalent time intervals) of the experiments
                   K_mean = c(sort(rnorm(5, 600, 350), decreasing = TRUE), # predictor variable (mean potassium of soil plots in treatment group)
                              sort(rnorm(5, 550, 325), decreasing = TRUE),
                              sort(rnorm(5, 650, 400), decreasing = TRUE)), 
                   biomass_effect = c( rnorm(5, 0.01, 0.3),  # response variable (biomass effect size)
                                       rnorm(5, -0.08, 0.4), 
                                       rnorm(5, 0.005, 0.2))) 

I have tried lme in the nlme package in R, with plant species as a random effect. However, this does not account for dependencies within the predictor variable.

install.packages("nlme")
library(nlme)
install.packages("car")
library(car) # for Anova

model.lme <- lme(biomass_effect ~ K_mean, random = ~1 | species, 
             correlation = corAR1(form = ~ round | species), data = data, method = "REML")

Anova(model.lme)

Another potential model moves the plant species from a random effect to a main effect with gls (in the nlme package).

model.gls <- gls(biomass_effect ~ K_mean + species, 
             correlation = corAR1(form = ~ round | species), data = data, method = "REML")

Anova(model.gls)

I am also interested in other potential ways to perform this analysis without summarizing the variables. I am concerned about losing variability associated with the response and predictor variables by using summary variables (effect size and mean) in the model.

$\endgroup$
  • $\begingroup$ I would not model this with an auto-correlation structure (difficult to do with three time points anyway). Instead I would include the soil plot and the rounds as random effects. $\endgroup$ – Roland Jun 15 '18 at 12:35
  • $\begingroup$ There are 5 time points (rounds) in this example, but that may still be difficult if you think 3 time points is difficult. Because I am using summary variables (mean K and effect size), the soil plots have been summarized and could not be used as a random effect in the model. $\endgroup$ – s.loftus Jun 15 '18 at 16:03
  • $\begingroup$ Well, don't use summary variables. $\endgroup$ – Roland Jun 15 '18 at 17:18
  • $\begingroup$ I'd like to not use summary variables, but am unsure how to represent the response without one, since I need to normalize the treatment response by the control response. To retain individual responses, would you recommend dividing each treatment soil plot response by the average control plot response? $\endgroup$ – s.loftus Jun 15 '18 at 19:23
  • $\begingroup$ Can't you include the control response as a predictor on the RHS instead of this normalization? $\endgroup$ – Roland Jun 15 '18 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.