2

I am trying to convert SAS to R for a mixed model of Repeated measures.The SAS code is as follows:

proc mixed data=pd method=reml ; 
by set; 
class id month arm;
model eff=base arm month arm*month/s;
repeated/subject=id type=un r;
lsmeans arm*month/pdiff cl;

In R, I have a data frame set up with month,id,arm as factors and eff,base as numeric. I need use the same method as above to generate the pvalues for each treatment arm (arm) and each month to see if their is a significant difference in the treatments during each given month.

My current attempt to replicate this in R is:

eff.lm<-lme(data = pd, fixed= eff~base + arm*month, 
    random= ~month|ID,
    correlation= corSymm(form=~month|ID))
eff.ls<-lsmeans(object= eff.lm, spec = ~arm|month)
comp<-pairs(eff.ls)
pval<-summary(pairs(eff.ls))$p.value

However, this results in the following error following the call to lme():

Error in if (length(uCov) != maxCov){  :
    missing value where TRUE/FALSE needed
In addition: Warning message:
In Ops.factor(unlist(covar),1):'-' not meaningful for factors

In response, I changed Month to Numeric and the model ran with no errors, but pairs function only generated pvalues at month 3.5 or the mean month, when I need pvalues for each individual month (1 to 6).

pd$month<-as.numeric(pd$month)
eff.lm<-lme(data = pd, fixed= eff~base + arm*month, 
    random= ~month|ID,
    correlation= corSymm(form=~month|ID))
eff.ls<-lsmeans(object= eff.lm, spec = ~arm|month)
comp<-pairs(eff.ls)
pval<-summary(comp)$p.value
print(comp)


> month = 3.5:
 contrast   estimate        SE   df t.ratio p.value
 1 - 2    -1.6817998 0.2415967 2015  -6.961  <.0001
 1 - 3    -1.2286784 0.2415994 2015  -5.086  <.0001
 1 - 4    -0.9165385 0.2416426 2015  -3.793  0.0029
 1 - 5    -0.7736447 0.2416409 2015  -3.202  0.0234
 1 - 6    -0.6709280 0.2416238 2015  -2.777  0.0809
 1 - 7    -0.1179885 0.2416134 2015  -0.488  0.9990
 2 - 3     0.4531214 0.2415986 2015   1.876  0.4969
 2 - 4     0.7652613 0.2416463 2015   3.167  0.0262
 2 - 5     0.9081551 0.2416446 2015   3.758  0.0033
 2 - 6     1.0108718 0.2416266 2015   4.184  0.0006
 2 - 7     1.5638113 0.2416156 2015   6.472  <.0001
 3 - 4     0.3121399 0.2416678 2015   1.292  0.8560
 3 - 5     0.4550337 0.2416657 2015   1.883  0.4919
 3 - 6     0.5577504 0.2416437 2015   2.308  0.2404
 3 - 7     1.1106899 0.2416296 2015   4.597  0.0001
 4 - 5     0.1428938 0.2415966 2015   0.591  0.9971
 4 - 6     0.2456105 0.2415991 2015   1.017  0.9503
 4 - 7     0.7985500 0.2416039 2015   3.305  0.0168
 5 - 6     0.1027167 0.2415987 2015   0.425  0.9995
 5 - 7     0.6556562 0.2416032 2015   2.714  0.0953
 6 - 7     0.5529395 0.2415979 2015   2.289  0.2499

P value adjustment: tukey method for comparing a family of 7 estimates

If anyone can offer suggestions on how to replicate the SAS code in R using with either the NLME package or the lme4 package, it would be appreciated.

Shayne03
  • 129
  • 4
  • 2
    Sidenote: the choices "lm" and "ls" are not optimal names for object as both of these represent functions for linear modelling and listing the names of objects in the environment respectively. This can lead to confusing error messages among other issues down the road. – lmo May 16 '17 at 14:00
  • 2
    any chance of a reproducible example? – Ben Bolker May 16 '17 at 14:27
  • also ... at what point does the error occur? only at the last step? – Ben Bolker May 16 '17 at 17:07
  • The error occurs following the call to lme() when month is represented as a factor. It does not occur when month is represented as numeric. The issue when month is numeric is that I cannot get the pairwise comparison at each individual month. It only outputs the pairwise comparison for month 3.5. – Shayne03 May 16 '17 at 17:35
  • @BenBolker What would be needed for a reproducible example? Sorry, I am new to using this site. – Shayne03 May 16 '17 at 17:38
  • http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – Ben Bolker May 16 '17 at 18:50
  • You can use `lsmeans(object= eff.lm, spec = ~arm|month, at = list(month=1:6))` as documented. But your model is not the same as your SAS model, which has `month` as a `class` variable., – Russ Lenth May 17 '17 at 02:53

0 Answers0