fitted values in mmrm

[davood1983]

Summary

Your brief description of the problem

# your reproducible example here

R session info

# R -e "utils::sessionInfo()" output goes here

OS / Environment

• OS: [e.g. Windows 10, Ubuntu 20.04, Centos 8]
• Docker Image [e.g. rocker/verse:4.1.0]

[danielinteractive]

Hi @davood1983 , can you please describe the question/problem?

[davood1983]

Hi, Thanks for getting back to me. I have a clinical experiment with 6 assessments. The first assessment is before applying the treatment so the means in the first assessment are equal. The model that I fit is mm2=mmrm(l.colony~ f.day+f.day:Treatment+ar1h(f.day|f.tu), data=d4) and then change the correlation matrix. This model imposes the equality of the means at assessment 1, however, when I check the fitted(mm2) this condition is not imposed. The other problem is that when I change the correlation matrix, the fitted values remain constant. Is mmrm compatible with the fitted(.) command? If not, what is your suggestion for getting the fitted values?

Best, Davood

[danielinteractive]

Hi Davood,

cool thanks for the initial information. Can you please provide a reproducible example, e.g. with fake data and the full code so I can run this and have a look?

The mmrm package supports the fitted method, see https://openpharma.github.io/mmrm/latest-tag/reference/mmrm_tmb_methods.html . Please also have a look at emmeans to obtain least square means, I am not sure if that is what you are looking for: https://openpharma.github.io/mmrm/latest-tag/articles/introduction.html?q=emmeans#support-for-emmeans

[davood1983]

library(mmrm)

# Set seed for reproducibility
set.seed(123)

# Number of levels
num_levels <- 4  # 3 treatment levels + 1 control

# Number of replicates for each level
num_replicates <- 5

# Number of assessments
num_assessments <- 6

# Mean values for each level at the first assessment
initial_means <- rep(10, num_levels)

# Simulate data
data <- data.frame(
  Subject = rep(1:(num_replicates * num_levels), each = num_assessments),
  Assessment = rep(1:num_assessments, times = num_replicates * num_levels),
  Treatment = rep(1:num_levels, each = num_replicates),
  Value = numeric(num_replicates * num_levels * num_assessments)
)

# Set initial values for each level
data$Value[data$Assessment == 1] <- rep(initial_means, each = num_replicates)

# Introduce variability to simulate treatment effect
treatment_effect <- c(0, 2, -1, 1)  # Adjust these values as needed
for (i in 2:num_assessments) {
  for (j in 1:num_levels) {
    data$Value[data$Assessment == i & data$Treatment == j] <- 
      data$Value[data$Assessment == (i - 1) & data$Treatment == j] + 
      rnorm(num_replicates, mean = treatment_effect[j], sd = 2)
  }
}

# Print the first few rows of the simulated dataset
head(data)
data$Assessment<-as.factor(data$Assessment)
data$Treatment<-as.factor(data$Treatment)
data$Subject<-as.factor(data$Subject)
n <- length((levels(data$Treatment)))

library(randomcoloR)
palette <- distinctColorPalette(n)

interaction.plot(
  x.factor = data$Assessment,
  trace.factor = data$Treatment,
  response = data$Value,
  col = palette,
  type = "b"
)
mm33 <- mmrm(Value~  Assessment+Assessment:Treatment+cs(Assessment|Subject), data=data)
summary(mm33)
data$c1 <- fitted(mm33)

d5 <- data[data$Assessment=="1",]

library(dplyr)
d5 |> 
  group_by(Treatment) |> 
  summarise(mean(c1))

interaction.plot(
  x.factor = data$Assessment,
  trace.factor = data$Treatment,
  response = c1,
  col = palette,
  type = "b"
)

I have a question, when the first level of the assessment is considered as the reference group, why there is an interaction between the reference group of the assessment and the treatment?

I think the reference group gets the zero value and the interactions should not be estimated in output. in my real data, these interactions are big numbers.

[danielinteractive]

Thanks @davood1983 , I just ran your example.

One comment: In the doubly nested loop I get 40 warnings:

There were 40 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In data$Value[data$Assessment == (i - 1) & data$Treatment ==  ... :
  longer object length is not a multiple of shorter object length

so it might well not do you want.

Also I needed to add some packages that were not included, I edited the code accordingly above. However c1 is not defined so the last statement does not pass.

If you would like a good answer from us, please ensure the example is really reproducible (run it in a fresh R session - does it work?) and please clearly describe the problem that you see or the question that you have. At the moment I still don't see what is the issue.

[danielinteractive]

Hi @davood1983 - if this is still an issue, please let me know, if I can run your example I am sure we can fix it!

[davood1983]

Thank you for your response. I used another package (nlme) even though it is an old package. I think that the problem is not related to the package itself. When we use time as a factor and remove the main effect of the treatment for baseline adjustment (response~Time(factor)+treatment(factor):Time(factor)) the interactions of the reference category of time and treatments remain in the model and are estimated. Based on the theory of linear models (since we do not have the sum(treatment=0)) this model is correct, but it is not the model that adjusts the model for the based line. However, when the time is a continuous variable and is scaled to get the zero value at baseline, none of these problems happen.

So, when time is a factor(discrete), we should add treatment to the model manually as dummies. Maybe it is worth mentioning it in your very helpful guideline for longitudinal data analysis.

In general, there is no problem in the package or model, it was a lack of knowledge about the terms in the model when time is a factor or continuous.

Again, thank you for your answer.

Best, Davood

On Thu, Feb 22, 2024 at 2:55 PM Daniel Sabanes Bove < @.***> wrote:

Hi @davood1983 https://github.com/davood1983 - if this is still an issue, please let me know, if I can run your example I am sure we can fix it!

— Reply to this email directly, view it on GitHub https://github.com/openpharma/mmrm/issues/421#issuecomment-1960166955, or unsubscribe https://github.com/notifications/unsubscribe-auth/AOZIS54OCTNAX5FAYNJBXILYU6PCNAVCNFSM6AAAAABCLHHKV6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNRQGE3DMOJVGU . You are receiving this because you were mentioned.Message ID: @.***>

[danielinteractive]

Sounds good, thanks @davood1983 !