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

I have the same question. ... View more

Is the correct interpretation of PROC GENMOD with a poisson distribution an odds ratio such as in logistic regression, or is the interpretation the same as an OLS Parameter Estimate? ... View more

For GLM: To test for homogeneity of variance across CLASS levels in a one-way ANOVA setting, use the MEANS statement with the HOVTEST option. There are a variety of tests available. I prefer the Brown-Forsythe method.   For mixed models (GLIMMIX being preferred); Include a RANDOM statement with the following syntax:   RANDOM _residual_/group = <class variable>; And add the following COVTEST statement:   COVTEST 'common variance' homogeneity This conducts a likelihood ratio test. One nice thing about this approach is that even for moderate sized datasets, you can ignore the testing part, and rely on convergence in distribution of the variance components if, in fact, there is no heteroscedasticity. You do this at the cost of degrees of freedom, so if you are in a situation where power of the test is at a premium, you can remove the statements if the test is non-significant.   For a situation where you have a main linear regression and you are looking for the influence of fixed categorical covariates, the GLIMMIX approach is quite handy.   SteveDenham     ... View more