trying to understand emmeans [General Sta­tis­tics]

posted by mittyri – Russia, 2019-11-02 02:26  – Posting: # 20745
Views: 1,921

Dear Zizou. dear All,

» the reason of the difference is that sequences were unbalanced (not equal number of subjects in each of the sequences).

looks like that reason is not the dominant reason.

» For unbalanced sequences there might be more questions ...

when the party comes to incomplete cases, everyone has more fun!

I am lazy man too, but some code for your pleasure

library(dplyr)
library(nlme)
library(replicateBE)
library(emmeans)
messupdataset<- function(dataset, treatment = 'R'){
  dataset$logPK <- ifelse(dataset$treatment==treatment, dataset$logPK*runif(1),dataset$logPK)
  return(dataset)
}

emmeanscomparison <- function(dataset, messuprefdata = F){
  dataset$logPK <- log(dataset$PK)
  if(messuprefdata==T){
    dataset <- messupdataset(dataset, "R")
  } 
 
  M=lm(logPK ~ sequence + subject + period + treatment, data = dataset)
  cat(paste0("\nLSMeans ratio by lm(): \nT/R = ", exp(coef(M)[["treatmentT"]])))
 
  newdat <- expand.grid(treatment = levels(dataset$treatment),
                        sequence = levels(dataset$sequence),
                        subject = levels(dataset$subject),
                        period = levels(dataset$period))
  preddata <- cbind(newdat, predict(M, newdat))
  preddatawoNA <- na.omit(left_join(preddata, dataset, by = c("treatment", "sequence", "subject", "period")))
  lsmeansbyhand <-
    preddatawoNA%>%
    group_by(treatment, sequence, period) %>%
    summarize(subjectmean = mean(`predict(M, newdat)`)) %>%
    group_by(treatment) %>%
    summarize(lsmean = exp(mean(subjectmean))) %>%
    mutate(lsmeansratio = lsmean[treatment=='T']/lsmean[treatment=='R']) %>%
    as.data.frame()
 
  cat(paste0("\nLSMeans by hand: \nT = ", lsmeansbyhand$lsmean[lsmeansbyhand$treatment=='T']))
  cat(paste0("\nLSMeans by hand: \nR = ", lsmeansbyhand$lsmean[lsmeansbyhand$treatment=='R']))
  cat(paste0("\nLSMeans ratio by hand: \nT/R = ", lsmeansbyhand$lsmeansratio[1], "\n"))
 
  emmeansdf <- data.frame(emmeans(M, 'treatment'))
  cat(paste0("\nLSMeans by emmeans: \nT = ", exp(emmeansdf$emmean[emmeansdf$treatment=='T'])))
  cat(paste0("\nLSMeans by emmeans: \nR = ", exp(emmeansdf$emmean[emmeansdf$treatment=='R'])))
  cat(paste0("\nLSMeans ratio by emmeans: \nT/R = ", exp(emmeansdf[emmeansdf$treatment=='T',]$emmean - emmeansdf[emmeansdf$treatment=='R',]$emmean), "\n"))
}

# balanced complete
emmeanscomparison(rds11, messuprefdata = F)
seed = 123
emmeanscomparison(rds11, messuprefdata = T)

# unbalanced complete
emmeanscomparison(rds16, messuprefdata = F)
seed = 123
emmeanscomparison(rds16, messuprefdata = T)

# balanced incompete
emmeanscomparison(rds26, messuprefdata = F)
seed = 123
emmeanscomparison(rds26, messuprefdata = T)

# unbalanced incompete
emmeanscomparison(rds14, messuprefdata = F)
seed = 123
emmeanscomparison(rds14, messuprefdata = T)


so I can get the same results for complete cases, but not for INcomplete. Please also take a look at the lsmeans for T when R is messed up for unbalanced complete and unbalanced incompltete. Looks like competeness has a predominant value
so my understanding of emmeans is also INcomplete and Unbalanced:-D

Kind regards,
Mittyri

Complete thread:

Activity
 Admin contact
20,801 posts in 4,354 threads, 1,446 registered users;
online 12 (0 registered, 12 guests [including 6 identified bots]).
Forum time: 04:45 CEST (Europe/Vienna)

Medical statistician: One who will not accept that Columbus discovered America…
because he said he was looking for India in the trial plan.    Stephen Senn

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5