R help - Feb 2015 - AR1 covariance structure for lme object and R/SAS differences in model output

Dear R users, 
We are working on a data set in which we have measured repeatedly a
physiological response variable (y) 
every 20 min for 12 h (time variable; 'x') in subjects ('id')
beloning to
one of five groups ('group'; 'A' to 'E'). Data are
located at:
https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0

We are interested to model if the response in y differences with time (i.e.
'x') for the two groups. Thus:
require(nlme)
m1<-lme(y~group*x+group*I(x^2),random=~x|id,data=data.df,na.action=na.omit)

But because data are collected repeatedly over short time intervals for each
subject, it seemed prudent to consider an autoregressive covariance
structure. Thus:
m2<-update(m1,~.,corr=corCAR1(form=~x|id))

AIC values indicate the latter (i.e. m2) as more appropriate:
  anova(m1,m2)
#   Model df      AIC      BIC       logLik        Test  L.Ratio     
p-value
#m1     1 19 2155.996 2260.767 -1058.9981                        
#m2     2 20 2021.944 2132.229  -990.9718 1 vs 2 136.0525  <.0001

Fixed effects and test statistics differ between models. A look at marginal
ANOVA tables suggest inference might differ somewhat between models:
  
anova.lme(m1,type="m")
#              numDF denDF  F-value p-value
#(Intercept)      1  1789 63384.80  <.0001
#group             4    45      1.29  0.2893
#x                   1  1789     0.05  0.8226
#I(x^2)            1  1789     4.02  0.0451
#group:x          4  1789     2.61  0.0341
#group:I(x^2)   4  1789     4.37  0.0016

anova.lme(m2,type="m")
#             numDF denDF  F-value p-value
#(Intercept)      1  1789 59395.79  <.0001
#group             4    45      1.33  0.2725
#x                   1  1789     0.04  0.8379
#I(x^2)            1  1789     2.28  0.1312
#group:x          4  1789     2.09  0.0802
#group:I(x^2)   4  1789     2.81  0.0244

Now, this is all well. But: my colleagues have been running the same data
set using PROC MIXED in SAS and come up with substantially different results
when comparing SAS default covariance structure (variance components) and
AR1. Specifically, there is virtually no change in either test statistics or
fitted values when using AR1 instead of Variance Components in SAS, which
fits the observation that AIC values (in SAS) indicate both covariance
structures fit data equally well. 

This is not very satisfactory to me, and I would be interesting to know what
is happening here. Realizing
this might not be the correct forum for this question, I would like to ask
you all if anyone would have any
input as to what is going on here, e.g. am I setting up my model
erroneously, etc.? 

N.b. I have no desire to replicate SAS results, but I would most certainly
be interested to know what could possibly explain  such a large discrepancy
between the two platforms. Any suggestions greatly welcomed.

(Data are located at:
https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0)

With all best wishes,
Andreas






--
View this message in context:
http://r.789695.n4.nabble.com/AR1-covariance-structure-for-lme-object-and-R-SAS-differences-in-model-output-tp4703103.html
Sent from the R help mailing list archive at Nabble.com.

Does anyone have code that uses a Kalman filter to impute time-series data? If
not, do you know of any software that can be used to impute time-series data?
Thank you,
John
John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric
Medicine
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing) 

Confidentiality Statement:
This email message, including any attachments, is for the sole use of the
intended recipient(s) and may contain confidential and privileged information.
Any unauthorized use, disclosure or distribution is prohibited. If you are not
the intended recipient, please contact the sender by reply email and destroy all
copies of the original message.

Hi,

searching for "kalman filter R" gives this paper that was published in
JSS. It may help.

@article{Tusell:2010:JSSOBK:v39i02,
  author =	"Fernando Tusell",
  title =	"Kalman Filtering in R",
  journal =	"Journal of Statistical Software",
  volume =	"39",
  number =	"2",
  pages =	"1--27",
  day =  	"1",
  month =	"3",
  year = 	"2011",
  CODEN =	"JSSOBK",
  ISSN = 	"1548-7660",
  bibdate =	"2010-08-17",
  URL =  	"http://www.jstatsoft.org/v39/i02",
  accepted =	"2010-08-17",
  acknowledgement = "",
  keywords =	"",
  submitted =	"2010-01-12",
}

15:23:00 -0500 "John Sorkin" <JSorkin at grecc.umaryland.edu>
wrote:
> Does anyone have code that uses a Kalman filter to impute time-series
> data? If not, do you know of any software that can be used to impute
> time-series data? Thank you, John
> John David Sorkin M.D., Ph.D.
> Professor of Medicine
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology and
> Geriatric Medicine Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing) 
> 
> Confidentiality Statement:
> This email message, including any attachments, is for ...{{dropped:28}}

> On 11 Feb 2015, at 16:57 , anord <andreas.nord at biol.lu.se> wrote:
> 
> Dear R users, 
> We are working on a data set in which we have measured repeatedly a
> physiological response variable (y) 
> every 20 min for 12 h (time variable; 'x') in subjects
('id') beloning to
> one of five groups ('group'; 'A' to 'E'). Data are
located at:
> https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0
> 
> We are interested to model if the response in y differences with time (i.e.
> 'x') for the two groups. Thus:
> require(nlme)
>
m1<-lme(y~group*x+group*I(x^2),random=~x|id,data=data.df,na.action=na.omit)
> 
> But because data are collected repeatedly over short time intervals for
each
> subject, it seemed prudent to consider an autoregressive covariance
> structure. Thus:
> m2<-update(m1,~.,corr=corCAR1(form=~x|id))
> 
> AIC values indicate the latter (i.e. m2) as more appropriate:
>  anova(m1,m2)
> #   Model df      AIC      BIC       logLik        Test  L.Ratio     
> p-value
> #m1     1 19 2155.996 2260.767 -1058.9981                        
> #m2     2 20 2021.944 2132.229  -990.9718 1 vs 2 136.0525  <.0001
> 
> Fixed effects and test statistics differ between models. A look at marginal
> ANOVA tables suggest inference might differ somewhat between models:
> 
> anova.lme(m1,type="m")
> #              numDF denDF  F-value p-value
> #(Intercept)      1  1789 63384.80  <.0001
> #group             4    45      1.29  0.2893
> #x                   1  1789     0.05  0.8226
> #I(x^2)            1  1789     4.02  0.0451
> #group:x          4  1789     2.61  0.0341
> #group:I(x^2)   4  1789     4.37  0.0016
> 
> anova.lme(m2,type="m")
> #             numDF denDF  F-value p-value
> #(Intercept)      1  1789 59395.79  <.0001
> #group             4    45      1.33  0.2725
> #x                   1  1789     0.04  0.8379
> #I(x^2)            1  1789     2.28  0.1312
> #group:x          4  1789     2.09  0.0802
> #group:I(x^2)   4  1789     2.81  0.0244
> 
> Now, this is all well. But: my colleagues have been running the same data
> set using PROC MIXED in SAS and come up with substantially different
results
> when comparing SAS default covariance structure (variance components) and
> AR1. Specifically, there is virtually no change in either test statistics
or
> fitted values when using AR1 instead of Variance Components in SAS, which
> fits the observation that AIC values (in SAS) indicate both covariance
> structures fit data equally well. 
> 
> This is not very satisfactory to me, and I would be interesting to know
what
> is happening here. Realizing
> this might not be the correct forum for this question, I would like to ask
> you all if anyone would have any
> input as to what is going on here, e.g. am I setting up my model
> erroneously, etc.? 
> 
> N.b. I have no desire to replicate SAS results, but I would most certainly
> be interested to know what could possibly explain  such a large discrepancy
> between the two platforms. Any suggestions greatly welcomed.
> 
> (Data are located at:
> https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0)
> 
Hmm, does SAS include a nugget effect perchance? At any rate, showing the SAS
output (or some of it) might make it easier for someone to help.

Also, R-sig-ME is a better choice for discussions of mixed effects models.


-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

anord <andreas.nord <at> biol.lu.se> writes:
> 

  [snip snip]
> We are working on a data set in which we have measured repeatedly a
> physiological response variable (y) 
> every 20 min for 12 h (time variable; 'x') in subjects
('id') beloning to
> one of five groups ('group'; 'A' to 'E'). Data are
located at:
> https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0
> 
> We are interested to model if the response in y differences 
> with time (i.e.
> 'x') for the two groups. Thus:
> require(nlme)
> m1<-lme(y~group*x+group*I(x^2),random=~x|id,
>   data=data.df,na.action=na.omit)
> 
> But because data are collected repeatedly over 
> short time intervals for each
> subject, it seemed prudent to consider an autoregressive covariance
> structure. Thus:
> m2<-update(m1,~.,corr=corCAR1(form=~x|id))
> 
> AIC values indicate the latter (i.e. m2) as more appropriate:
>   anova(m1,m2)
> #   Model df      AIC      BIC       logLik        Test  L.Ratio     
> p-value
> #m1     1 19 2155.996 2260.767 -1058.9981                        
> #m2     2 20 2021.944 2132.229  -990.9718 1 vs 2 136.0525  <.0001
> 
> Fixed effects and test statistics differ between models.
>  A look at marginal
> ANOVA tables suggest inference might differ somewhat between models:
> 
> anova.lme(m1,type="m")
> #              numDF denDF  F-value p-value
> #(Intercept)      1  1789 63384.80  <.0001
> #group             4    45      1.29  0.2893
> #x                   1  1789     0.05  0.8226
> #I(x^2)            1  1789     4.02  0.0451
> #group:x          4  1789     2.61  0.0341
> #group:I(x^2)   4  1789     4.37  0.0016
> 
> anova.lme(m2,type="m")
> #             numDF denDF  F-value p-value
> #(Intercept)      1  1789 59395.79  <.0001
> #group             4    45      1.33  0.2725
> #x                   1  1789     0.04  0.8379
> #I(x^2)            1  1789     2.28  0.1312
> #group:x          4  1789     2.09  0.0802
> #group:I(x^2)   4  1789     2.81  0.0244
> 
> Now, this is all well. But: my colleagues have been running 
> the same data
> set using PROC MIXED in SAS and come up with 
> substantially different results
> when comparing SAS default covariance structure (variance components) and
> AR1. Specifically, there is virtually no change 
> in either test statistics or
> fitted values when using AR1 instead of Variance Components in SAS, which
> fits the observation that AIC values (in SAS) indicate both covariance
> structures fit data equally well. 
> 
> This is not very satisfactory to me, and I would be 
> interesting to know what
> is happening here. Realizing
> this might not be the correct forum for this question, I would like to ask
> you all if anyone would have any
> input as to what is going on here, e.g. am I setting up my model
> erroneously, etc.? 
> 
> N.b. I have no desire to replicate SAS results, but I would most certainly
> be interested to know what could possibly explain  
> such a large discrepancy
> between the two platforms. Any suggestions greatly welcomed.
> 
> (Data are located at:
> https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0)
  Could you repost this on r-sig-mixed-models at r-project.org ?

It would be useful to see some of the comparisons (fixed effects,
RE variance-covariance, correlation parameter) between SAS and
R; are they actually fitting the same model?  (e.g., does SAS
allow for covariance between the slope and intercept random effects?)
Maybe they're getting the same estimates but computing df/p-values
in different ways?

  I thought I would try this with orthogonal polynomials in case
some of the fits were unstable ...

data.df <- read.csv2("ar1_data.csv")
library("nlme")
m1 <- lme(y~group*x+group*I(x^2),random=~x|id,
          data=data.df,na.action=na.omit,method="ML")
## use method="ML" so we can compare orthog. and non-orthog.
polynomials
m1B <- update(m1,.~group*poly(x,2))
m2 <- update(m1,corr=corCAR1(form=~x|id))
m2B <- update(m1B,corr=corCAR1(form=~x|id))
AIC(m1,m1B,m2,m2B)

Thanks for your comments. I hard disk meltdown has slowed me down somewhat,
but I am now back online. 

I have reposted to the ME-list, and requested SAS output from my colleagues.
I will update the thread again as soon I know more.

Thanks,
Andreas



--
View this message in context:
http://r.789695.n4.nabble.com/Autoregressive-covariance-structure-for-lme-object-and-R-SAS-differences-in-model-output-tp4703103p4703396.html
Sent from the R help mailing list archive at Nabble.com.