openpharma/mmrm

fitted values in mmrm

Closed this issue · 8 comments

davood1983 commented

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 commented

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

davood1983 commented

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 commented

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 commented 
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 commented

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 commented

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 commented
danielinteractive commented

Sounds good, thanks @davood1983 !