1

model suggested by FDA

I am struggling on how to fit the model for replicated crossover design using REML method. The suggested model by FDA is as above and can someone help on how to code it into R coding ? This is my coding and I wonder if it is right or wrong ?

samplePK.lmer <- lmer(ykir2~1+Treatment:Sequence:Replication+
    (1|Subject:Sequence:Treatment), data=samplePK, REML=TRUE)
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
fat
  • 21
  • 2
  • PS is this from here? https://juliahub.com/docs/ReplicateBE/RNURj/1.0.10/details/ It would be interesting to run comparisons of this package with `MixedModels.jl` (in Julia) and `lme4` (in R). – Ben Bolker Dec 21 '20 at 17:05
  • @BenBolker I have installed the julia package in r but unfortunately there are some errors. Can you please help me to find a way to fix it ? – fat Dec 23 '20 at 03:35
  • library("JuliaCall") Warning message: package ‘JuliaCall’ was built under R version 3.6.3 > julia_setup() Julia version 1.5.3 at location C:\Users\User\AppData\Local\Programs\JULIA1~1.3\bin will be used. Error in inDL(x, as.logical(local), as.logical(now), ...) : unable to load shared object 'C:/Users/User/AppData/Local/Programs/JULIA1~1.3/bin/libjulia.dll': LoadLibrary failure: %1 is not a valid Win32 application. Error in juliacall_initialize(.julia$dll_file) : C:\Users\User\AppData\Local\Programs\JULIA1~1.3\bin\libjulia.dll - %1 is not a valid Win32 application. – fat Dec 23 '20 at 03:58
  • I'm sorry; the Julia/R interface is frequently broken, I have a hard time keeping up with it myself. – Ben Bolker Dec 23 '20 at 17:33

1 Answers1

0

I would say the formula should be

response ~ trt + trt:seq:rep + (trt|subj:seq)

The key difference from your specification is that (trt|subj:seq) is fitting something like a randomized-block or random-slopes model, where we are allowing for the variation of the treatment effect across subject/sequence combinations.

There are a few issues here that I ran into/noticed when trying to fit this model to simulated data:

  • if we are fitting this with "modern" mixed-model approaches, there will be some parameters aliased between the treatment effect (trt) and the trt:seq:rep term. In a method-of-moments approach this doesn't matter so much (because you never explicitly estimate the parameters), but it leads to complaints about rank-deficient models (which are ignorable if you know what you're doing ...).
  • it seems wrong that the random effect delta_{ij} is given as having a mean of {mu_R,mu_T}; this is redundant with the fixed effect mu_k

Obviously I could have something wrong or misunderstood something about the original model specification.

I might suggest that you try follow-up questions on the r-sig-mixed-models@r-project.org mailing list, where there is a wide readership with broad expertise (i.e., more expertise on the subject of mixed models specifically than here on StackOverflow)

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • So I've changed the fitted model and yes it shows the complaint about the rank-deficient models. Is it okay if I received the message and can I just continue with the model ? – fat Dec 22 '20 at 04:53
  • I think so. It would be really nice to do some cross-checking across packages/platforms (SAS, ReplicateBE.jl, etc.) to be sure ... – Ben Bolker Dec 22 '20 at 16:44