在R中重新编码SAS PROC混合模型-Java 学习之路

我正在尝试将混合模型从SAS重新编码为R，而我正处于艰难时期，因为SAS不是我曾经使用过的编码语言 . 我一直试图用三个文件来帮助我：

但是，我仍然遇到适合该模型的问题 .

我正在使用的数据集来自 rmcorr 包

library(rmcorr)
library(dplyr)
data <- bland1995
data <- as.data.frame(data %>%
            group_by(Subject) %>%
            mutate(Test = seq_along(Subject)))

head(data)
  Subject   pH PacO2
1       1 6.68  3.97
2       1 6.53  4.12
3       1 6.43  4.09
4       1 6.33  3.97
5       2 6.85  5.27
6       2 7.06  5.37

SAS混合的目的是计算pH和PacO2之间的部分测量相关性，同时控制受试者的重复测量 . SAS代码如下：

proc mixed method = ml covtest;
class patient mvar replicate;
model response = mvar / solution ddfm = kr;
random mvar / type = un subject = patient v vcorr;
repeated mvar replicate/ type = un@ar(1) subject = patient r rcorr;
run;

在上面的代码中，patient = subject和mvar表示数据集，其格式与我当前拥有的数据集不同 . 我不确定这是否与SAS有关？当我尝试使用SAS使用的格式的数据对模型进行编码时，我仍然没有达到相同的结果 . SAS的数据采用长格式，如下所示：

library(reshape)
library(dplyr)

data <- as.data.frame(data %>%
                group_by(Subject) %>%
                mutate(Test = seq_along(Subject)))

data.new <- melt(data, id = c("Subject", "Test"), measure.vars = c("pH", "PacO2"))

data.new <- data.new[order(data.new$Subject, data.new$Test), ]

head(data.new)
   Subject Test variable value Test.f
1        1    1       pH  6.68      1
48       1    1    PacO2  3.97      1
2        1    2       pH  6.53      2
49       1    2    PacO2  4.12      2
3        1    3       pH  6.43      3
50       1    3    PacO2  4.09      3

我试图在R中重新编码这个模型的方式如下：

library(nlme)

m1 <- lme(pH ~ PacO2, random = ~Test|Subject, data = data, correlation = corAR1(form = ~Test|Subject))

summary(m1)

Linear mixed-effects model fit by REML
 Data: data 
        AIC       BIC   logLik
  -66.36773 -53.72109 40.18386

Random effects:
 Formula: ~Test | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev     Corr  
(Intercept) 0.26422171 (Intr)
Test        0.04809647 0.14  
Residual    0.09860705       

Correlation Structure: AR(1)
 Formula: ~Test | Subject 
 Parameter estimate(s):
      Phi 
0.7791259 
Fixed effects: pH ~ PacO2 
                Value  Std.Error DF  t-value p-value
(Intercept)  7.555977 0.13154582 38 57.43989       0
PacO2       -0.079537 0.01641498 38 -4.84536       0
 Correlation: 
      (Intr)
PacO2 -0.635

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.17707008 -0.40505338 -0.04571545  0.25761717  1.15654722 

Number of Observations: 47
Number of Groups: 8

SAS模型报道pH和PacO2之间的部分测量相关性为0.0464 . 我似乎无法理解这个数字是如何从模型中获得的（假设我的模型以相同的方式编码 - 我不确定它是什么） . 我试图将方差 - 协方差矩阵看作是理解这种相关性的一种方式，这是SAS模型用来得出答案的方法：

vcov(m1)
VarCorr(m1)
getVarCov(m1, type = "marginal", individual = 2)
getVarCov(m1, type = "random.effects", individual = 2)

但这似乎都不正确 . 有什么我想念的吗？我还可以提供我用来尝试在长格式数据集上创建模型的代码，如果有帮助的话？

在R中重新编码SAS PROC混合模型

相关问题