PROC MIXED - covariance structures for equally and unequally spaced time series

yaofengluan

请问SAS MIXED中对适用于等距和不等距时间序列的covariance structures 分别有哪些？

以下的见解对吗？

等距：type = cs ; type = toep ; type = ar(1)

不等距： type = vc; type = sp(pow)(c-list)

yaofengluan

我找了很多材料，但是，对于不等距的尤其模糊。请知道的同学和老师指教。

rtist

proc mixed的文档里面不是有两个表么，都列的很清楚了。

yaofengluan

回复第3楼的 rtist：

我又读了一遍那两个表及相关说明，但是还是没能理解，不确定到底那些可以归到等距，那些可以归到不等距。我认为Spatial Covariance Structures的这些可以归到不等距吧，但是像 type = vc 跟 type = cs 能不能用到不等距里面呢？？？请赐教，不甚感激。

rtist

I don't see why they couldn't be used, especially when there are very few time pionts. If you indeed have an infinitely long time series, probably they are not used so often.

yaofengluan

Thanks for your answer.

There are six time points in all.

This is not so many, so that I can use “type = cs" for testing unequally spaced time series. Is is right?

rtist

This is not a time series. Type=cs is a "possible" choice, not necessarily a good choice. "Usually", you can rewrite it as a random effect model to get the same structure. As of which structure to use, probably one needs extra background information (which is rarely available) or use some data dependent measures, like AIC or BIC (these are practical, but I don't see any formal justification unless there are a large number of independent blocks), to decide.

yaofengluan

回复第7楼的 Rtist：谢谢！我也是在别人的间接帮助下（那个老师请别的老师帮我），一步一步认识SAS的。对于MIXED model 尤其不甚了解。我用了 type = cs,而且与我期望的结果更靠近。但是，我不知道是不是对的。我介绍下我的试验，并将别人帮我写的程序给你看下。

我从田里采集昆虫，每15天左右采集一次，遇到雨天顺延。数据运用重复测量进行分析。

<br />
options ls=80 ps=55 nocenter pageno=1 mlogic mprint symbolgen;</p>
<p>%let path = %sysget(SAS_EXECFILEPATH);<br />
%global path;<br />
footnote "&path";</p>
<p>%macro names(name= , number= );<br />
  %do i = 1 %to &number;<br />
	  &name&i<br />
	%end;<br />
%mend names;</p>
<p>Data Fengluan;<br />
  input Landid Trt Blk N1 N2 N3 N4 N5 N6;<br />
	one = 1;<br />
  Time=1; TimeJJ = 178; N=N1; LogN=log10(N1); output;<br />
  Time=2; TimeJJ = 193; N=N2; LogN=log10(N2); output;<br />
  Time=3; TimeJJ = 208; N=N3; LogN=log10(N3); output;<br />
  Time=4; TimeJJ = 223; N=N4; LogN=log10(N4); output;<br />
  Time=5; TimeJJ = 238; N=N5; LogN=log10(N5); output;<br />
  Time=6; TimeJJ = 253; N=N6; LogN=log10(N6); output;<br />
	keep LandId Trt Blk one Time TimeJJ N LogN;<br />
  datalines;<br />
1  1  1  7  40  154  205  323  104<br />
2  1  2  3  15  305  176  119  79<br />
3  1  3  2  19  278  118  69  86<br />
4  1  4  6  36  109  123  71  86<br />
5  2  1  1  13  193  145  97  160<br />
6  2  2  4  17  218  159  20  35<br />
7  2  3  1  6  253  91  23  53<br />
8  2  4  2  15  171  94  30  69<br />
9  3  1  8  35  251  228  248  152<br />
10  3  2  3  16  365  193  107  60<br />
11  3  3  2  20  252  135  115  106<br />
12  3  4  7  18  334  131  35  74<br />
;</p>
<p>proc format;<br />
  value Trtshow 1='CK' 2='Co/R' 3='S/R';<br />
run;</p>
<p>proc format;<br />
  value Blkshow 1='A' 2='B' 3='C' 4='D';<br />
run;</p>
<p>proc format;<br />
  value Timeshow 1='06-27-09' 2='07-12-09' 3='07-27-09'<br />
                 4='08-11-09' 5='08-26-09' 6='09-10-09';<br />
run;</p>
<p>proc sort data=Fengluan;<br />
  by  Trt Blk  Landid Time;<br />
run;</p>
<p>PROC PRINT DATA=Fengluan;<br />
  Title 'DATA Fengluan';<br />
  format Trt Trtshow.;<br />
  format Blk Blkshow.;<br />
  format Time Timeshow.;<br />
run;</p>
<p>* Create 5 orthogonal polynomial contrasts to partition the effect of time;</p>
<p>proc iml;<br />
  TEMPS = {178, 193, 208, 223, 238, 253};<br />
  OT = orpol(TEMPS,5);<br />
  OTP = t(OT[,2:6]);<br />
  LABT = 'T1':'T5';<br />
  create TO from OTP [rowname=LABT];<br />
  append from OTP [rowname=LABT];<br />
  close TO;<br />
quit;<br />
data CTO;<br />
  set TO;<br />
  file 'C:\Users\Jarron\Desktop\CT.TXT';<br />
  put '  contrast  ''' LABT ''' Time ' / @8 (COL1 - COL6) (12.8) ';';<br />
run;</p>
<p>%macro mix(type);<br />
Proc mixed DATA=Fengluan;<br />
  title "Unequally spaced repeated measures: type=&type";<br />
  Class Trt Blk  Landid Time;<br />
  Model LogN= Blk Trt Time Trt*Time / outp=res DDFM=KR;<br />
  Repeated Time / sub=landid type=&type RCI;<br />
  %include 'C:\Users\Jarron\Desktop\CT.TXT';<br />
	contrast 'Trt 1 vs 2+3' Trt -2  1  1;<br />
	contrast 'Trt 2 vs 3' Trt  0 -1  1;<br />
	lsmeans Trt / cl adjust=simulate;<br />
	ods output InvCholR=RCI;<br />
	ods listing exclude InvCholR;<br />
run;<br />
proc iml;<br />
  use RCI;<br />
	read all var{%names(name=Col,number=6)} into RCI;<br />
	close RCI;<br />
	use res;<br />
	read all var{resid} into resid;<br />
  close res;</p>
<p>	* de-correlate residuals (more or less);</p>
<p>	res_cortd = (I(12)@RCI)*resid;</p>
<p>	create res2 from res_cortd [colname={'res_cortd'}];<br />
	append from res_cortd;<br />
quit;</p>
<p>data res_plus;<br />
  merge res res2;<br />
run;</p>
<p>proc univariate data=res_plus plot normal;<br />
  %if %bquote(&type) ~= VC %then %str(var resid res_cortd;);<br />
  %else %str(var resid;);<br />
run;</p>
<p>proc plot data=res_plus;<br />
  plot resid*(pred Trt Time Blk);<br />
  %if %bquote(&type) ~= VC %then %str(plot res_cortd*(pred Trt Time Blk););<br />
quit;<br />
%mend mix;</p>
<p>* %mix(un);  * Unreliable,  too many var-cov parms;<br />
%mix(sp(pow)(timeJJ one)); * slightly better than VC;<br />
%mix(VC);<br />
%mix(CS); *Fengluan added;<br />

</p>

yaofengluan

最后一行是我加上去的：%mix(CS); *Fengluan added;

yaofengluan

这边的是等距离的时间，我2008年的数据是不等距离的。

试验设计是随机区组设计，三个处理，四次重复。重复数不多，老师建议当做固定因子

rtist

Results from the codes seem indicating some heterogeneity in variance across the time points.

Depending the goal of your analysis, this may or may not be troublesome.

The CS results indicated negative correlation, which needs careful thinking about its interpretability.

The DDFM=KR option implemented in SAS is problematic for nonlinear covariance structures; see Kenward & Roger (2009, CSDA).

It'll be interesting to see the results from type=un hlps as well.

Two more interesting papers to read are Skene & Kenward (2010, StatMed).

If all you need is a test, I'll suggest trying some permutation procedures, as you sample size is too small to justify most asymptotic tests.

yaofengluan

回复第11楼的 Rtist：谢谢，我不是统计学专业的，我学的是昆虫，尽管统计非常重要。但是，这些内容对于目前的我还是太深奥了。

CS 算出的结果是我期望的，但是，这并不一定符合客观。如果我用GLM来算的话，结果也是我所期望的。

我用我期望的结果写成了文章，准备投稿，不知道会遇到啥情况。呵呵，谢谢你的热心帮助。

顺便说下，你的英语好好啊，羡慕。我用英文写的话，很费劲，而且不知道你会不会看的懂，就用中文了。

投稿后，遇到新问题，我直接给你留言啊