Dear @SteveDenham, Thank you again for your response. 1. Stage 1 or 2 My code is for stage 2, stage 1 appeared to work. I've posted both stage 1 and stage 2 of my code below with comments for clarification. Maybe something notable for stage 1: if I used proc mixed without the model statement: singular=1e-8 I will occasionally get the "Infinite Likelihood" error. I never get the error otherwise or in glimmix. I don't think this is due to duplicates, even though that's normally the case. I got the idea to use singular=1e-8 from a thread post by @jiltao. https://communities.sas.com/t5/SAS-Programming/Infinite-Likelihood-error-Proc-Mixed/m-p/803817#M316523 Other than scouring the data by hand, I used tabulate to make sure I only have one observation per treatment combination per year. proc tabulate data=dry_full;
class year amend cut rep;
table (amend*cut*rep)*(N), year ALL; run; 2. REML You're right, of course. The code I pasted isn't using REML, but I have been using REML up until recently. I've entered the "just try anything" phase of the analysis, and changing the method was the latest thing. 3. The last of the 42 must be the residual variance and is not included in the code (yes/no ?) The 42nd covariance parameter is included in the code and set to 1 (hold=42). It isn't obvious from the code I posted originally because I had imported parmsdata from a previous run. It should be more clear below. The large blocks of parms values have 20 values each (20 treatment combinations). The first block of 20 are for the Shukla stability intercepts, and the second block of 20 are for the Shukla stability slopes. There is one additional parms value before the blocks (random year) and one after (residual). 4. Is there anything in the Macholdt paper that explains the values that they plug into the stage 1 parms statement? Not to my knowledge, and I don't think it's always necessary to set parms values in stage 1. It looks to me as though they wanted to guarantee non-negative AR(1) values since there isn't a "nobounds" option anywhere in their code. This wasn't a problem in my data. Here's another appendix screenshot, from a different paper. Note that they're also using glimmix in stage 2, but they aren't estimating Shukla slopes. Another appendix 5. Results from splitting the data for stage 2 into two parts: estimates/treatments and the covariance matrix The results are identical when the data are split, SAS does not seem to be confused by having the same dataset passed to "ldata" and "data". 6. Results from specifying parms values based on the hacky regression estimates I re-ran the data with very specific parms estimates that were based exactly on (1) a linear model of the per-decade Shukla estimates and (2) a quadratic model of the same estimates. The linear model results in: ERROR: Values given in PARMS statement are not feasible. The quadratic model has insufficient memory, so I'll re-run it this weekend (we have SAS studio with an 8gb limit... I can get 16gb over the weekend if I ask nicely). 7. Full code with comments and explanations as to parm bounds. Stage 1 /* I do not use the parms statement at all in stage 1.*/ /* This is a split plot experiment (cut main, amend sub)*/ /* I used SP(POW) for the autocorrelation in the first two random statements*/ /* because we have "missing" (fallow) years (unequal time steps)*/ /* I used csh for the autocorrelation in the third random statement because*/ /* I needed to account for heteroscedasticity more so than an AR1 structure (IMO).*/ /*SP(POW) doesn't seem to work nicely with heterogeneous variances.*/ /* CSH seemed to be a good compromise since other structures (UN, ANTE) */ /* resulted in insufficient memory errors*/ Proc glimmix Data=DRY_full ;
ods output covparms=cp2 lsmeans=dryyield_full_csh;
Class rep year amend cut;
Model tr1_yield=cut|amend|year /ddfm=kr ;
random year/subject=rep type=SP(POW)(year_num);
random year/subject=rep*cut type=SP(POW)(year_num);
random year/subject=rep*cut*amend type=csh;
Lsmeans cut*amend*year/cov;
run;
/* The data are modified to get new variables "t", "st", and conform to the lin(1) naming conventions*/
Data dryyield_full_csh_changed;
Set dryyield_full_csh;
Row=_N_;
Rename cov1-cov572=col1-col572;
Parm=1;
t=year-1990;
st=sqrt(t);
Run; Stage 2 Proc glimmix Data=dryyield_full_csh_changed plots=pearsonpanel(conditional marginal) method=rspl;
Class year amend cut cult ;
Model estimate=amend|cut|t cult /s ddfm=kr cl noint;
nloptions maxiter=1000 GCONV=0;
COVTEST /cl;
Random year; /*Results in 1 covariance est*/
Random int/sub=year group=amend*cut; /*20 covariance est*/
Random st/sub=year group=amend*cut; /*20 covariance est*/
random _residual_/subject=int type=lin(1) ldata=dry_full; /*1 covariance est*/
Parms /*This is the "year" starting value*/ (.01)
/*These are the Shukla intercept starting values*/ /*My guesses presented here came from the "hacky" regression*/ /*The two 0.00000000 values may be an issue, */ /*the regression estimated them as being negative. */ /*Removing the lower bound did not help for them.*/ /*The quadratic version may actually help if so (all were positive)*/ /*There are 20 intercept values, 20 slope values*/
(0.06804396)(0.00000000)(0.03140187)(0.11450479)
(0.11665396)(0.16136292)(0.00823989)(0.00894313)
(0.02657896)(0.02634313)(0.24657283)(0.05338596)
(0.03234812)(0.03043496)(0.06711858)(0.11232925)
(0.01128913)(0.00000000)(0.03359046)(0.03478354)
/*These are the Shukla slope starting values*/
(-0.00332221)(00.00006651)(-0.00180324)(-0.00668694)
(-0.00559049)(-0.00787263)(-0.00045049)(-0.00053659)
(-0.00149999)(-0.00078679)(-0.01335398)(-0.00303279)
(-0.00172644)(-0.00167246)(-0.00325155)(-0.00501563)
(-0.00021841)(00.00032859)(-0.00184381)(-0.00143606)
/*This is the covariance parameter associated with the ldata covariance matrix*/ /*It is held at unity as per the papers I've read on the this 2-stage method*/ /*i.e. hold=42, it is the 42nd value*/
(1)/
hold=42
nobound
lowerb= /*The variances have a lower bound of zero*/ 0,
0,0,0,0,
0,0,0,0,
0,0,0,0,
0,0,0,0,
0,0,0,0, /*The slopes for the Shukla stability estimates can be negative*/
.,0,.,.,
.,.,.,.,
.,.,.,.,
.,.,.,.,
.,0,.,.,
/*This value is held at 1, so a lowerbound probably doesn't matter...*/
0;
lsmeans amend*cut;
Run; This above code results in "parms not feasible". If we comment out the "nobound" option and use (0.01) for all parms except lin(1), we get: NOBOUND commented out of the code If we then import these covariance estimates (parmsdata=cp2) and include NOBOUNDS we get: Generic parms values, NOBOUND active I'll try the quadratic version over the weekend when I have more memory. Thanks very much for taking a look at this, it's been driving me crazy...
... View more
