these look good to me. My only suggestion would be on the multiple years without a repeated statement. Consider rewriting the RANDOM statement with block(site) as a subject. You would then have a random intercept and random slope model:   proc mixed ; class Year Site Block Treatment; model biomass = Year|Site|Treatment; random intercept Treatment/subject=Block(Site) type=un; run; I added type=un to handle any potential correlation between the intercept and Treatment.  This one does bother me a bit. I don't have an issue with a fixed effect being modeled as an R-side effect (repeated measures approach), but I think you will run into convergence and/or G matrix issues with Treatment as both random and fixed.   SteveDenham     ... View more

We all tend to see different things, but it appears that the absolute value of the residuals increase with increasing predicted values. You might consider that to be evidence of heteroskedasticity, but it might also be the result of model misspecification, improper scaling of design variables, collinearity (if this is a multiple regression) or not recognizing that whatever process generates the data is one where scale is a function of location, such as a gamma distribution.   SteveDenham ... View more

you might want to search Mr. Google for repeated AB testing, and see how others have analyzed their results.   SteveDenham ... View more

I was able to run your model without errors, but there are some results that make me wonder if your model could be over specified. Here is the code I used:   PROC GLIMMIX DATA = test1/*(where=(stratclus<999))*/ METHOD=QUAD(FASTQUAD QPOINTS=20) plots=all; CLASS RAHHIDPN STRATCLUS GENDER(REF='1') RACE_ETH(ref='1') MARITAL(ref='1'); PARMS (0.01) (1.72) (0.5); MODEL MOBILITY = BSAGE_C BSAGE_C*BSAGE_C GENDER GENDER*BSAGE_C RACE_ETH RACE_ETH*BSAGE_C RACE_ETH*GENDER TIME TIME*TIME TIME*BSAGE_C TIME*GENDER TIME*RACE_ETH/solution obsweight=HCNSWGTR ;/**LEVEL 1 VARIABLES**/ RANDOM INTERCEPT /subject=STRATCLUS type=VC ; RANDOM INTERCEPT /subject=RAHHIDPN(STRATCLUS) type=VC ; NLOPTIONS GCONV=0.00001 tech=congra; COVTEST/WALD; ODS OUTPUT parameterestimates=FIXED1 covparms=RAND1; RUN; The two changes were to reduce the number of quadrature points, and to use the PARMS statement to provide some better starting values for the random effects.  Regarding my thought of over-specification, the F values for five of the terms in the model went to infinity. These were GENDER, RACE_ETH*GENDER, TIME*TIME, TIME*RSAGE_C and TIME*GENDER. You may want to look to see if there is something unusual with the strata and the number of observations available in the various cells.   SteveDenham   ... View more

Thank you for your responses. My data is actually monthly, but I have dropped the ‘Month’ column. I have tried the ADF test as a secondary test when validating the model. I have tried proc arima but don’t know how to set the AIC/BIC option for lag length. This is my code: proc reg data=df; model dep = ind1 ind2 ind3 ; OUTPUT OUT=model_output_c predicted= pred residual=resid; run; proc arima data=model_output_c; identify var=resid STATIONARITY=(ADF) ; run; ... View more

Further, how can a duration be binomially distributed? Unless you are defining an incidence as a duration, such as behavior not seen --> duration=0, behavior seen -> 1.   But your issue is with the first analysis, and the problem is that the log of 0 is negative infinity, so folks say you can't analyze the results when you have a zero and a log link. There are alternatives, and you are part way there by picking a multiplicative overdispersion. Consider adding a small (say 1e-6) to every observation, and fitting a gamma distribution (although you can still fit a Poisson, let's just move to the continuous version rather than staying with the discrete version). That should address your model. It is harder to compare means except on the log scale. The ILINK option gives a ratio for the exponentiated difference of two means on the log scale, not the difference on the measured scale of two exponentiated estimates. I don't know if the %NLest macro would be appropriate here, but it certainly is appropriate for the Poisson distribution example. Perhaps @StatDave could chime in here on that approach.   SteveDenham ... View more

The syntax would be the same as for the above design, so I will repeat that code, and then insert your factors.   /* previous poster */ PROC GLIMMIX data=data; class site origin block cut tree; model recovery = site|origin|cut/ddfm=kr; random block block*site block*origin block*cut block*site*origin block*site*cut block*origin*cut block*origin*site*cut;; /* change variable names to organic, system, nitrogen */ PROC GLIMMIX data=data; class block organic system nitrogen; model recovery = organic|system|nitrogen/ddfm=kr2; random block block*organic block*system block*nitrogen block*organic*system block*organic*nitrogen block*organic*system*nitrogen; Again, I want to point out that you probably do not have enough data to estimate the three-way and four-way effects, and they probably should be considered residual error in any case. To improve the ability to converge, you may then want this for the RANDOM statement:   random intercept organic system nitrogen/subject=block; SteveDenham     ... View more

Thanks for your explanations. I am actually t modeling raw BP values over time (BP1 (baseline) and BP2(first follow-up). And Var1 is typo. ... View more

Thanks Mr Steve Denham, you are Olympian to me I checked the book, BUGS and your solution and I love it THANKS     ... View more

Use an ODS OUTPUT statement to get the ESTIMATE value into a dataset, and then postprocess that dataset, comparing the result of the ESTIMATE statement to -1. This might be easier if you use LSMESTIMATE as you don"t need all of the terms that refer to anything except the lsmeans for visit*baseline.   SteveDenham ... View more

You may need to simplify your model, as I see 3 variables that are functions of one another:   Chg = observed value - base Base Strata = median function determined by Base,   so the induced collinearity between base and strata may well be driving the similarity of the estimates between the strata groups. What happens if you remove strata from the model, and (perhaps) post-process the predicted values obtained from an OUTP= option in your MODEL statement into quintiles, terciles or two groups? Does that address your research question?   SteveDenham ... View more

Like @SteveDenham mentioned, you can create a new variable indicating which type the response is, turn the data into a long format and fit a model accordingly. Please note that there is only one MODEL statement allowed in PROC MIXED, and only one response variable in the MODEL statement to be specified. Here is an example to show you this approach: subject  height weight gender 1              165         56          F 2             180         76          M .... data two;  set one;  response=height;  type='Height'; output; response=weight; type='Weight'; output; drop height weight; run; Then use Response as your dependent variable in the MODEL statement, and include TYPE as one of your effect in the MODEL statement. In your case, the values for TYPE would be your original response variables. Hope this helps, Jill ... View more

The GLIMMIX documentation goes through several caveats regarding the ordering of the categories for a nominal multinomial analysis.  In the end, I think their approach would look like this when translated to your situation:   proc glimmix data=lab.merge INITGLM; class occ_reduced region(ref='SK') ; model occ_reduced(order= freq ref='LOW') = diffH diffM diffL region / solution DIST=MULTINOMIAL link=glogit; where labour=1; weight weight; run; Quoting the documentation: " In generalized logit models (for multinomial data with unordered categories), one response category is chosen as the reference category in the formulation of the generalized logits. By default, the linear predictor in the reference category is set to 0, and the reference category corresponds to the entry in the "Response Profile" table with the highest Ordered Value. You can affect the assignment of Ordered Values with the DESCENDING and ORDER= options in the MODEL statement. You can choose a different reference category with the REF= option. The choice of the reference category for generalized logit models affects the results. It is sometimes recommended that you choose the category with the highest frequency as the reference (see, for example, Brown and Prescott 1999, p. 160). You can achieve this with the GLIMMIX procedure by combining the ORDER= and REF= options..."   I don't know if that will help or not, but my experience with generalized logit models is that critical clues to its function in GLIMMIX aren't stressed enough in the examples.   SteveDenham ... View more

I would caution you about any inferences regarding drug or treatment. The design I see completely confounds drug*treatment with pen, and for a lot of folks that would mean that pen is the experimental unit, and animals are observational units measured on each of the pens.  You might be able to use the following, but I provide no guarantees. I moved to the spatial power covariance structure because I believe the residuals will be correlated with those closer in time more highly correlated. The autoregressive error structures (including ANTE and ARMA) treat the time points as being equally spaced in time. I also fit pen as an effect, rather than trt or drug. You can get LSMEAN-like estimates for the various levels of each by using LSMESTIMATE statements to manipulate the pen lsmeans. I added the /e option to the LSMEANS statement to facilitate construction of the LSMESTIMATES:   proc mixed data= XX; class animal pen day; model BW = pen|day initial_BW ddfm=kr; random intercept / subject = animal; repeated days / subject = animal type = sp(pow)(t); /* t is a continuous version of 'day' */ lsmeans pen|day /e adjust = tukey; /* Add a section of LSMESTIMATE code here to calculate trt and drug means and differences */ run; SteveDenham   ... View more

If you have 20G observations, and have DEN F coming out as Inf, think about what that means - you are out at the point where F tests are indistinguishable from chi squared tests, and the t distribution doesn't differ from the normal until you get out to several decimal places. So, rather than trying DDFM= basically anything, try using the CHISQ option, and use those results rather than the F tests.   If you ran the analysis in GLIMMIX, you could use ddfm=none.   SteveDenham ... View more