[R-sig-ME] Converting SAS code into R

[Viechtbauer, Wolfgang (NP)]

Dear Nicholas,

The model specified via 'proc mixed' is actually not a mixed-effects / multilevel model at all. The default 'type' for 'repeated' is VC, which is the same as in a standard regression type model (i.e., it assumes independence between repeated observations within subjects and assumes a single variance error variance component). So you are in essence just fitting:

proc glm data=rmanova4; class randomization_arm cancer_type site wk;
 model chgpf = randomization_arm cancer_type site wk;
run;

which is identical in R to:

res <- lm(chgpf ~ Randomization_Arm + Cancer_Type + site + wk, data=rmanova.data)
summary(res)

(assuming that the the predictors are all coded as factors in R).

I suspect that you actually wanted to fit a different model in SAS to begin with (why else use proc mixed?), but you would have to tell us first what the actual model is that you wanted to fit in SAS before we can tell you how to fit the same model in R.

Best,
Wolfgang

-- 
Wolfgang Viechtbauer, PhD, Statistician | Department of Psychiatry and    
Neuropsychology | Maastricht University | PO Box 616 (VIJV1) | 6200 MD    
Maastricht, The Netherlands | +31(43)3884170 | https://www.wvbauer.com    

>-----Original Message-----
>From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-project.org] On
>Behalf Of Nicholas Mitsakakis via R-sig-mixed-models
>Sent: Thursday, 05 January, 2023 17:45
>To: r-sig-mixed-models using r-project.org
>Subject: [R-sig-ME] Converting SAS code into R
>
>I have been trying to convert some PROC MIXED SAS code into R, but without
>success. The code is:
>
>proc mixed data=rmanova4;class randomization_arm cancer_type site wk;
>model chgpf=randomization_arm cancer_type site wk;
>repeated  / subject=study_id;
>contrast '12 vs 4' randomization_arm 1 -1;
>lsmeans randomization_arm / cl pdiff alpha=0.05;
>run;quit;
>
>I have tried something like
>
>lme(chgpf ~ Randomization_Arm  + Cancer_Type + site + wk, data=rmanova.data,
>            random = ~ 1  | Study_ID,
>            correlation = corSymm(form = ~ 1 | Study_ID),
>            na.action=na.exclude, method = "REML")
>
>but I am getting different estimate values.
>
>Perhaps I am misunderstanding something basic. Any comment/suggestion would
>be greatly appreciated.
>
>I am adding here the output. Part of the output from the SAS code is below:
>
>Least Squares Means
>Effect  Randomization_Arm   Estimate    Standard Error  DF  t Value Pr
>> |t|    Alpha   Lower   Upper
>Randomization_Arm   12 weekly BTA   -4.5441 1.3163  222 -3.45   0.0007
> 0.05    -7.1382 -1.9501
>Randomization_Arm   4 weekly BTA    -6.4224 1.3143  222 -4.89   <.0001
> 0.05    -9.0126 -3.8322
>
>Differences of Least Squares Means
>Effect  Randomization_Arm   _Randomization_Arm  Estimate    Standard
>Error  DF  t Value Pr > |t|    Alpha   Lower   Upper
>Randomization_Arm   12 weekly BTA   4 weekly BTA    1.8783  1.4774
>222 1.27    0.2049  0.05    -1.0332 4.7898
>
>The output from the R code is below:
>
>Linear mixed-effects model fit by REML
>  Data: rmanova.data
>       AIC     BIC    logLik
>  6526.315 6586.65 -3250.157
>
>Random effects:
> Formula: ~1 | Study_ID
>        (Intercept) Residual
>StdDev:    16.51816 12.95417
>
>Correlation Structure: General
> Formula: ~1 | Study_ID
> Parameter estimate(s):
> Correlation:
>  1      2      3     2  0.222              3 -0.159  0.225       4
>-0.421 -0.042  0.083
>Fixed effects:  chgpf ~ Randomization_Arm + Cancer_Type + site + wk
>                                  Value Std.Error  DF    t-value
>p-value(Intercept)                    5.739240 2.8987216 541
>1.9799209  0.0482
>Randomization_Arm4 weekly BTA -1.174704 2.3915873 225 -0.4911817  0.6238
>Cancer_TypeProsta             -4.459715 2.4711025 225 -1.8047469  0.0725
>site                          -1.917902 0.9709655 225 -1.9752530  0.0495
>wk                            -1.570707 0.5115174 541 -3.0706809  0.0022
> Correlation:
>                              (Intr) R_A4wB Cnc_TP site
>Randomization_Arm4 weekly BTA -0.440
>Cancer_TypeProsta             -0.314  0.043
>site                          -0.598  0.003 -0.064
>wk                            -0.421 -0.004  0.032  0.003
>
>Standardized Within-Group Residuals:
>        Min          Q1         Med          Q3         Max
>-4.99530346 -0.36507852  0.08308708  0.45937864  3.12730244
>
>Number of Observations: 771
>Number of Groups: 229
>
>--
>Nicholas Mitsakakis, MSc, PhD, P.Stat.
>
>Senior Biostatistician and Associate Scientist
>Children's Hospital of Eastern Ontario Research Institute
>Adjunct Lecturer, Dalla Lana School of Public Health, University of Toronto