用R的nlme::lme复现SAS PROC MIXED含随机效应与相关结构的模型

TL;DR 我正尝试将SAS的PROC MIXED模型用R的nlme::lme实现，但遇到困难。

我希望将Hamlet等人（2003）文献中的SAS代码复现为R代码。作者使用的SAS代码基于PROC MIXED：

proc mixed;
class persnum vtype replicate;
model response = vtype / solution ddfm=kr;
random vtype / type=un subject=persnum g gcorr v vcorr;
repeated vtype / type=un subject=replicate(persnum) r rcorr;
run;

我知道solution ddfm=kr指定了计算固定效应分母自由度的Kenward-Roger方法，该方法在R中不可用。忽略这一点，我尝试用nlme::lme()（甚至geepack::geeglm()，尽管我知道这会替换随机效应）复现模型。我的R代码如下：

model <-  nlme::lme(
    fixed = response ~ vtype
    ,random = list( ID = pdSymm( ~ vtype) )
    ,correlation = corSymm( form = ~ vtype | persnum/replicate)
    ,method = "ML"
    ,data = phPACO_long
  )

模型运行报错：

Error in if (length(uCov) != maxCov) { : missing value where TRUE/FALSE needed

我在nlme的Git仓库中找不到uCov或maxCov，因此不知道如何定位并修复该错误。有人了解这个错误的来源及解决方法吗？

我看到TueBimet仅通过nlme::lme()的random参数，复现了带有REPEATED语句但无RANDOM语句的PROC MIXED代码。但我没有SAS，无法验证该方法的有效性。此外，我要转换的SAS代码中RANDOM和REPEATED语句的分组不同，因此我认为该方法不适用于我的情况。

长格式数据的结构如下：

structure(list(persnum= structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L), levels = c("1", "2", "3", "4", "5", 
"6", "7", "8"), class = "factor"), replicate= structure(c(1L, 1L, 
2L, 2L, 3L, 3L, 4L, 4L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 1L, 1L, 
2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 
1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 1L, 1L, 2L, 2L, 3L, 3L, 
4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 1L, 1L, 2L, 2L, 3L, 3L, 
4L, 4L, 5L, 5L, 6L, 6L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 
3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L), levels = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9"), class = "ordered"), vtype= structure(c(2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), levels = c("pH", "P_aCO2"), class = "factor"), response= c(6.68, 3.97, 6.53, 4.12, 
6.43, 4.09, 6.33, 3.97, 6.85, 5.27, 7.06, 5.37, 7.13, 5.41, 7.17, 
5.44, 7.4, 5.67, 7.42, 3.64, 7.41, 4.32, 7.37, 4.73, 7.34, 4.96, 
7.35, 5.04, 7.28, 5.22, 7.3, 4.82, 7.34, 5.07, 7.36, 5.67, 7.33, 
5.1, 7.29, 5.53, 7.3, 4.75, 7.35, 5.51, 7.35, 4.28, 7.3, 4.44, 
7.3, 4.32, 7.37, 3.23, 7.27, 4.46, 7.28, 4.72, 7.32, 4.75, 7.32, 
4.99, 7.38, 4.78, 7.3, 4.73, 7.29, 5.12, 7.33, 4.93, 7.31, 5.03, 
7.33, 4.93, 6.86, 6.85, 6.94, 6.44, 6.92, 6.52, 7.19, 5.28, 7.29, 
4.56, 7.21, 4.34, 7.25, 4.32, 7.2, 4.41, 7.19, 3.69, 6.77, 6.09, 
6.82, 5.58), persnum.replicate= structure(c(1L, 1L, 9L, 9L, 17L, 17L, 
25L, 25L, 2L, 2L, 10L, 10L, 18L, 18L, 26L, 26L, 3L, 3L, 11L, 
11L, 19L, 19L, 27L, 27L, 32L, 32L, 37L, 37L, 41L, 41L, 44L, 44L, 
47L, 47L, 4L, 4L, 12L, 12L, 20L, 20L, 28L, 28L, 33L, 33L, 5L, 
5L, 13L, 13L, 21L, 21L, 29L, 29L, 34L, 34L, 38L, 38L, 42L, 42L, 
45L, 45L, 6L, 6L, 14L, 14L, 22L, 22L, 30L, 30L, 35L, 35L, 39L, 
39L, 7L, 7L, 15L, 15L, 23L, 23L, 8L, 8L, 16L, 16L, 24L, 24L, 
31L, 31L, 36L, 36L, 40L, 40L, 43L, 43L, 46L, 46L), levels = c("1.1", 
"2.1", "3.1", "4.1", "5.1", "6.1", "7.1", "8.1", "1.2", "2.2", 
"3.2", "4.2", "5.2", "6.2", "7.2", "8.2", "1.3", "2.3", "3.3", 
"4.3", "5.3", "6.3", "7.3", "8.3", "1.4", "2.4", "3.4", "4.4", 
"5.4", "6.4", "8.4", "3.5", "4.5", "5.5", "6.5", "8.5", "3.6", 
"5.6", "6.6", "8.6", "3.7", "5.7", "8.7", "3.8", "5.8", "8.8", 
"3.9"), class = c("ordered", "factor"))), row.names = c(NA, -94L
), class = c("tbl_df", "tbl", "data.frame"))

附言：如果在correlation参数中排除mvar，模型可以拟合，但得到的协方差矩阵与Hamlett等人的结果不匹配（这符合预期，因为它假设了不同的相关结构）。

内容的提问来源于stack exchange，提问作者Ciarán D. McInerney