Перекодирование SAS PROC MIXED MODEL в R

Question

Я пытаюсь перекодировать смешанную модель из SAS в R, и мне трудно, так как SAS не является языком кодирования, с которым я когда-либо работал. Я пытался использовать три документа, чтобы помочь мне:

Резюме смешанных моделей в SAS и R

lmer для смешанных пользователей SAS PROC

Определение линейных смешанных моделей в статистических пакетах

Тем не менее, у меня все еще есть проблемы с подгонкой модели.

Набор данных, который я использую, взят из 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;

В приведенном выше коде пациент = субъект, а 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)

Но все это не кажется правильным. Есть что-то, чего мне не хватает? Я также могу предоставить код, который использовал, чтобы попытаться создать модель поверх набора данных в длинном формате, если это могло бы помочь?