Using repeated measures change/trend as a predictor variable

#1
I have repeated measures of happiness for a sample of participants, and a single measure of satisfaction for each of the participants.

I want to predict satisfaction from the repeated measures of happiness. To do so, I want to create a new variable, called happiness.change, which measures the degree of change/trend in happiness from the first measurement to the last, for each participant (for example, a negative happiness.change if there is a decrease in happiness over time). Then I want to predict satisfaction from happiness.change.

Below (using R) is an excerpt from my data (a sample of 9 participants):

Code:
ids <- c(rep(seq(1:5), each = 2), rep(6:9, each = 5))
time <- c(rep(1:2, 5), rep(1:5, 4))
happiness <- c(0.80,0.00,0.75,0.00,0.80,0.00,2.75,2.50,0.40,0.20,
               3.80,0.40,0.00,0.20,3.40,3.00,2.60,3.40,3.80,0.00,
               3.60,3.60,0.20,0.40,1.00,0.40,0.20,1.20,1.20,0.00)
satisfaction <- c(6,6,2,2,3,3,2,2,2,2,5,5,5,5,5,7,7,7,7,7,1,1,1,1,1,2,2,2,2,2)
data <- as.data.frame(matrix(c(ids, time, happiness, satisfaction),
                         nrow = 30,
                         ncol = 4,
                         dimnames = list(c(),c("id", "time",
                                               "happiness", "satisfaction"))))
print(data)

#    id time happiness satisfaction
# 1   1    1      0.80            6
# 2   1    2      0.00            6
# 3   2    1      0.75            2
# 4   2    2      0.00            2
# 5   3    1      0.80            3
# 6   3    2      0.00            3
# 7   4    1      2.75            2
# 8   4    2      2.50            2
# 9   5    1      0.40            2
# 10  5    2      0.20            2
# 11  6    1      3.80            5
# 12  6    2      0.40            5
# 13  6    3      0.00            5
# 14  6    4      0.20            5
# 15  6    5      3.40            5
# 16  7    1      3.00            7
# 17  7    2      2.60            7
# 18  7    3      3.40            7
# 19  7    4      3.80            7
# 20  7    5      0.00            7
# 21  8    1      3.60            1
# 22  8    2      3.60            1
# 23  8    3      0.20            1
# 24  8    4      0.40            1
# 25  8    5      1.00            1
# 26  9    1      0.40            2
# 27  9    2      0.20            2
# 28  9    3      1.20            2
# 29  9    4      1.20            2
# 30  9    5      0.00            2
To create happiness.change I was advised to use the coefficients produced by either of these equations:

Code:
require(lme4)
model1 <- lmer(happiness ~ time + (time | id), data = data)

require(nlme)
model2 <- lme(happiness ~ time, random = ~1 + time | id, data = data)
For example, running coef(model1) produces the following coefficients (column time):

Code:
# $id
#   (Intercept)        time
# 1   0.9936158 -0.05991770
# 2   0.9739595 -0.05674569
# 3   0.9936158 -0.05991770
# 4   2.5539086 -0.31170766
# 5   0.8998610 -0.04478815
# 6   2.0446747 -0.22953078
# 7   3.2906564 -0.43059926
# 8   2.7120530 -0.33722798
# 9   1.0027053 -0.06138450
# 
# attr(,"class")
# [1] "coef.mer"
And then connecting between satisfaction and the coefficients:

Code:
coefs1 <- as.data.frame(unlist(coef(model1))[10:18])
satisfactionData <- reshape(data,
                            direction = "wide",
                            idvar = "id",
                            timevar = "time")[c(1,3)]
newData <- cbind(satisfactionData, coefs1)
colnames(newData) <- c("id", "satisfaction", "happiness.change")
print(newData)

#    id satisfaction happiness.change
# 1   1            6      -0.05991770
# 3   2            2      -0.05674569
# 5   3            3      -0.05991770
# 7   4            2      -0.31170766
# 9   5            2      -0.04478815
# 11  6            5      -0.22953078
# 16  7            7      -0.43059926
# 21  8            1      -0.33722798
# 26  9            2      -0.06138450
----------------------------------------

I have several questions:

  1. Running model2 generates the following error:

    Code:
    Error in lme.formula(happiness ~ time, random = ~1 + time | id, data = data) : 
      nlminb problem, convergence error code = 1
      message = iteration limit reached without convergence (10)
    I don't remember where, but I read that running lme(happiness ~ time, random = ~1 + time | id, control = list(opt = "optim"), data = data) bypasses this error, and indeed it does. But what exactly does it do? Do the coefficients produced by using this model still represent the change in happiness over time?

  2. Running model2 (with control = list(opt = "optim")) produces slightly (very slightly) different coefficients than model1; why? What is the difference between the models?

  3. What would be the suitable method for testing the relationship between happiness.change and satisfaction? I tried cor.test(newData$satisfaction, newData$happiness.change), which produced the following results, but I'm not sure how to interpret them (they are significant in the complete dataset):

    Code:
        Pearson's product-moment correlation
    
    data:  newData$satisfaction and newData$happiness.change
    t = -0.72735, df = 7, p-value = 0.4906
    alternative hypothesis: true correlation is not equal to 0
    95 percent confidence interval:
     -0.7901064  0.4843018
    sample estimates:
           cor 
    -0.2650786
  4. As you can see, model1 produces negative coefficients for all the participants (model2 as well). Even for participants whose change over time is purely positive (for example, a participant with 2 happiness measures: time1 = 0 and time2 = 0.8; no such example in the sample above, but it is in my data).

    This may be a problem, because I want to be able to distinguish between participants whose change in happiness is positive (happiness increases over time) and participants whose change in happiness is negative (happiness decreases over time); and then see whether there's a difference between these participants in their relationship between happiness.change and satisfaction. So my question is this: Is it statistically "legit" to divide my participants beforehand into groups based on their "raw" change in happiness (for example, I would label a participant with time1 = 0 and time2 = 0.8 as having a positive change), and then model the coefficient for these groups separately?

    However, creating sub-groups this way may be difficult with participants with more than 2 measures of happiness, which brings me to my final question:

  5. If I understand correctly, model1 and model2 assume there is a linear change over time. However, I checked the whole sample (424 participants), and a cubic model (happiness ~ time) actually explains more of the variance in happiness. So I suppose a cubic change over time is more appropriate. How can I create a new variable reflecting such a change in happiness over time?

I realize this is quite long and these are a lot of questions, and so any help will be greatly appreciated. If anyone can answer even one of these questions, I will be very grateful.

Thanks!