SAS proc glm随机效应模型的对比转化为R

Question

我对任何错误表示歉意；我直到最近才开始学习SAS。我得到了这个SAS代码(下面的代码是一个reprex，而不是确切的代码)，它使用proc glm来建立随机效应模型。SAS代码没有使用color，而是使用对比和idnumber间接映射到color。

我想知道如何在R中复制这一点。几次使用lme4进行随机效果的尝试和使用MASS::ginv的对比都失败了，所以我可能需要使用我不熟悉的包。

我还想知道red-blue和red-blue2之间的区别，以及为什么输出是不同的。谢谢你的帮助。

data df1;
input idnumber color $ value1;
datalines;
1001 red 189
1002 red 145
1003 red 210
1004 red 194
1005 red 127
1006 red 189
1007 blue 145
1008 red 210
1009 red 194
1010 red 127.
;

proc glm data=df1;
class idnumber;
    model value1=idnumber/noint solution clparm;
    contrast 'red vs. blue' idnumber 1 1 1 1 1 1 -9 1 1 1;
    estimate 'red-blue' idnumber 1 1 1 1 1 1 -9 1 1 1/ divisor=10;
    estimate 'red-blue2' idnumber .111 .111 .111 .111 .111 .111 -.999 .111 .111 .111;
run;

下面是几次复制尝试。

idnumber <- c(1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010)
color <- c('red', 'red', 'red', 'red', 'red', 'red', 'blue', 'red', 'red', 'red')
value1 <- c(189, 145, 210, 194, 127, 189, 145, 210, 194, 127)
df1 <- data.frame(idnumber, color, value1)

library(lme4)
library(MASS)
library(tidyverse)
options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))

# attempt 1
mod1 <- lme4::lmer(value1 ~ (1|idnumber), data = df1) # error
#> Error: number of levels of each grouping factor must be < number of observations (problems: idnumber)

# attempt 2
mod2 <- lme4::lmer(value1 ~ (1|color), data = df1) # singular
#> boundary (singular) fit: see help('isSingular')
summary(mod2)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: value1 ~ (1 | color)
#>    Data: df1
#> 
#> REML criterion at convergence: 90.9
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -1.3847 -0.8429  0.4816  0.6321  1.1138 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  color    (Intercept)    0      0.00   
#>  Residual             1104     33.22   
#> Number of obs: 10, groups:  color, 2
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)   173.00      10.51   16.47
#> optimizer (nloptwrap) convergence code: 0 (OK)
#> boundary (singular) fit: see help('isSingular')

# attempt 3
mat1 <- rbind(c(-0.5, 0.5))
cMat1 <- MASS::ginv(mat1)
mod3 <- lm(value1 ~ color, data = df1, contrasts = list(color = cMat1))
summary(mod3)
#> 
#> Call:
#> lm(formula = value1 ~ color, data = df1, contrasts = list(color = cMat1))
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -49.11 -23.33  12.89  17.89  33.89 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   160.56      17.74   9.052 1.78e-05 ***
#> color1         15.56      17.74   0.877    0.406    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 33.65 on 8 degrees of freedom
#> Multiple R-squared:  0.08771,    Adjusted R-squared:  -0.02633 
#> F-statistic: 0.7691 on 1 and 8 DF,  p-value: 0.4061

# attempt 4
con <- c(.1, .1, .1, .1, .1, .1, -.9, .1, .1, .1)
mod4 <- lm(value1 ~ idnumber, data = df1, contrasts = list(idnumber = con)) # error, but unsure how to fix
#> Error in `contrasts<-`(`*tmp*`, value = contrasts.arg[[nn]]): contrasts apply only to factors
Created on 2022-02-08 by the reprex package (v2.0.1)

Answer 1

回答这部分是可以回答的:这两种不同的估计是怎么回事。

估计语句包括系数列表。这些值乘以值，然后进行求和--给出结果。他们不一样的原因是他们不一样..。第一个是(除后) 0.1 / -0.9，第二个是0.111 (1/9)/ -0.999，实际上与第一个除数为9而不是10的除数相同。因此，数学是不同的。

我也不确定您的reprex，使用idnumber作为类变量没有任何意义.您似乎更可能使用color作为类变量。这有可能只是糟糕的SAS代码吗？我不是GLM专家，但尝试使用GLM时使用GLM却很奇怪，因为分类变量是ID号(假设它是唯一的ID )。