fitted values in mmrm
Closed this issue · 8 comments
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]
Hi @davood1983 , can you please describe the question/problem?
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
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
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.
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.
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!
Sounds good, thanks @davood1983 !