lme() in R - Linear mixed model for repeated-measures design. Fit independent slopes to determine associations within levels of the repeated measure?

Question

I did an experiment in which 30 people listened to 3 songs ("A", "B", and "C") in a random order. For each participant and song, I got the mean "body temperature" and the mean "heart rate" during the listening task. This is, I have a repeated-measures design with 1 factor (song) and 3 levels.

Here is some dummy data in long format following the structure of my real data:

res_table <- tibble(
  ID                   = rep(seq(1,30,1),3),
  SONGS                = rep(c("A","B","C"),each=30),
  BodyTemperature      = sample(35:65, 30*3, replace = TRUE),
  HeartRate            = sample(50:100, 30*3, replace = TRUE),
  YearsMusicalTraining = rep(sample(0:7,30, replace = TRUE),3)
)
res_table$ID <- factor(res_table$ID, labels = c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"))
res_table$SONGS <- factor(res_table$SONGS, levels = c("A","B","C"))

I am interested in knowing the degree of association between "body temperature" and "heart rate" across songs. Specifically, I want to know how well can the measures of "heart rate" predict "body temperature" and whether this association changes depending on the song the participants were listening to.

I tried to make a mixed model / multilevel linear model with R by using lme() by following Field's book "Discovering Statistics Using R":

model.01 <- lme(BodyTemperature ~ HeartRate, random = ~1|ID/SONGS, data = res_table, method = "ML")
summary(model.01)

but without much success: the summary of the ouput does not give me statistics of the associations between "body temperature" and "heart rate" for each song level. This is the output (I did not include AIC and BIC here but it is also part of the output):

Linear mixed-effects model fit by maximum likelihood
  Data: res_table 

Random effects:
 Formula: ~1 | ID
         (Intercept)
StdDev: 0.0005227154

 Formula: ~1 | SONGS %in% ID
        (Intercept)   Residual
StdDev:    9.189513 0.04996696

Fixed effects:  BodyTemperature ~ HeartRate 
 Correlation: 
          (Intr)
HeartRate -0.979

Standardized Within-Group Residuals:
          Min            Q1           Med            Q3           Max 
-0.0094852307 -0.0046700058  0.0001771823  0.0048029654  0.0084200939 

Number of Observations: 90
Number of Groups: 
           ID SONGS %in% ID 
           30            90 

Second part of the output

I also thought about building this model:

model.02 <- lme(BodyTemperature ~ HeartRate*SONGS, random = ~1|ID/SONGS, data = res_table, method = "ML")
summary(model.02)

but I cannot fully comprehend what the output of the summary means, even if I create contrasts beforehand to "put labels" on the created dummy variables:

AvB <- c(1,-1,0) # comparison of A versus B
CvB  <- c(0,1,-1) # comparison of C versus B
contrasts(res_table$SONGS)<-cbind(AvB,CvB)

This is the output (I did not include AIC and BIC here but it is also part of the output):

Linear mixed-effects model fit by maximum likelihood
  Data: res_table 

Random effects:
 Formula: ~1 | ID
         (Intercept)
StdDev: 0.0005119628

 Formula: ~1 | SONGS %in% ID
        (Intercept)   Residual
StdDev:    8.914116 0.04965589

Fixed effects:  BodyTemperature ~ HeartRate * SONGS 
 Correlation: 
                   (Intr) HertRt SONGSA SONGSC HR:SONGSA
HeartRate          -0.980                               
SONGSAvB            0.124 -0.103                        
SONGSCvB            0.030 -0.031  0.552                 
HeartRate:SONGSAvB -0.104  0.081 -0.981 -0.534          
HeartRate:SONGSCvB -0.031  0.033 -0.527 -0.979  0.527   

Standardized Within-Group Residuals:
          Min            Q1           Med            Q3           Max 
-0.0100366260 -0.0048459431  0.0002849439  0.0045537959  0.0103488878 

Number of Observations: 90
Number of Groups: 
           ID SONGS %in% ID 
           30            90 

Second part of the output

I do not understand the output of the summary of model.02 because I still only have the interaction between HeartRate and SONGSAvB, or between HeartRate and SONGSCvB (this is, between HeartRate in general and the two contrasts). Also, should I be looking at the values of the correlation to assess the association between the variables of interest across levels? Or at the coefficients of the model?

So, to be more specific, I do not understand: (1) how to build a model to look for the degree of association between the two variables of interest for each song level (am I not building the correct model? what is the correct formula?), and (2) which coefficient determines the degree of association between variables of interest (would the Mean Square of the model be a good indicator of this? Should I look at the values of the output where it says "Correlation"?)

I suspect I should create a regression line with a different intercept and slope for each song level but within the repeated-measures design. I have the intuition that building a linear model for each song level independently would be incorrect because of the multiple comparisons problem and because there is a within-subject variance that should be accounted for by the model. However, I do not where to start and I cannot find anything similar on the internet.

I hope my code and questions are clear enough eventhough I am quite new with R and statistics.

Thank you in advance,

R.