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

问题描述

我的目标是对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))
}

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

非常感谢您。

马可

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）