如何获得线性混合效应模型系数的最大似然估计和支持间隔？

如何解决如何获得线性混合效应模型系数的最大似然估计和支持间隔？

我的目标是对R中的线性混合效应模型使用似然范式（Royall，1997）。

因此，我想获得线性混合效应模型系数的最大似然估计（MLE）和支持间隔（SI）（1 / 6.8、1 / 8、1 / 20、1 / 32）。

这是我的数据集的一瞥。

Rows: 310
Columns: 8
$ name             <int> 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,1,2...
$ date             <date> 2019-07-17,2019-07-17,2019-0...
$ playing_position <chr> "Prop","Prop","Second Row","Hooker","Back Row",...
$ time             <dbl> 1,1...
$ therapist        <int> 1,2...
$ total_distance_m <dbl> 233.8404,257.5557,209.8772,217.8274,225.8800,253.0800,254.604...
$ player_load      <dbl> 26.63188,34.72477,30.13677,25.55523,25.60363,30.42476,33.0126...
$ meterage_min     <dbl> 93.53615,103.02229,83.95090,87.13096,90.35198,101.23198,101.8...

我在nlme中创建的模型是 my_model <- lme(player_load ~ time + therapist,random = ~ 1 + therapist | name,data = df)

我找到了ProfileLikelihood软件包，使用该模型alternative_model <- lme(player_load ~ time + therapist,random = ~ 1 | name,data = df)，我用它确定了时间系数和治疗师的MLE和SI。

mylike <- profilelike.lme(player_load ~ therapist,random = ~ 1 | id,data = df1,profile.theta = "time",method = "ML",subject = "name",lo.theta = -1,hi.theta = 6,length = 500,round = 2)
profilelike.plot(theta=mylike$theta,profile.lik.norm=mylike$profile.lik.norm,round=4)

但是，它不允许我使用随机斜率（即therapist）和计算截距的MLE和SI。

这是ProfileLikelihood::profilelike.lme()函数的代码，为了考虑随机斜率和截距系数的计算，我想对其进行一些修改。

function (formula,data,subject,random,correlation = NULL,profile.theta,lo.theta,hi.theta,length = 300,round = 2,subset = NULL,weights = NULL,...) 
{
  if (!is.null(subset)) {
    stop("Warning message: 'subset' should not be provided")
  }
  if (!is.null(weights)) {
    stop("Warning message: 'weights' should not be provided")
  }
  m <- model.frame(formula,data)
  X <- model.matrix(formula,m)
  y <- model.response(m)
  theta.off <- data[,names(data) == profile.theta]
  id <- data[,names(data) == subject]
  if (!is.numeric(theta.off)) {
    stop("Warning message: 'profile.theta' must be a numeric variable")
  }
  if ((length(theta.off) != length(y) | length(theta.off) != 
       length(X[,1]) | length(y) != length(X[,1]))) {
    cat("Warning message: remove missing data \n")
  }
  if ((is.null(lo.theta) | is.null(hi.theta))) {
    cat("Warning message: provide lo.theta and hi.theta \n")
    fit <- lm(y ~ -1 + X + theta.off,na.action = na.fail)
    mle <- summary(fit)$coefficient["theta.off",1]
    se <- summary(fit)$coefficient["theta.off",2]
    lo.theta <- round(mle - 4 * se,round)
    hi.theta <- round(mle + 4 * se,round)
  }
  theta <- seq(from = lo.theta,to = hi.theta,length = length)
  log.lik <- rep(NA,length)
  for (i in 1:length) {
    pi <- theta[i]
    y.off <- y - pi * theta.off
    fit <- lme(y.off ~ -1 + X,random = random,correlation = correlation,na.action = na.fail)
    log.lik[i] <- logLik(fit)
  }
  theta <- theta[is.na(log.lik) != 1]
  log.lik <- log.lik[is.na(log.lik) != 1]
  profile.lik <- exp(log.lik)
  mm <- max(log.lik,na.rm = TRUE)
  log.norm.lik <- log.lik - mm
  profile.lik.norm <- exp(log.norm.lik)
  return(list(theta = theta,profile.lik = profile.lik,profile.lik.norm = profile.lik.norm))
}

如果有人有任何想法，我将不胜感激，将不胜感激。

非常感谢您。

马可