Hello, I am performing logistic regression in PROC GLIMMIX (binomial distribution, logit link, SAS version 9.4) using categorical predictors and am coming across a strange phenomenon: If the response variable of certain level of a categorical predictor comprioses mostly or all zero values (my data are represented by 0s ands 1s -I am modelling the probability of occurrence of 1s), when I perform lsmean estimates of the probability of occurrence of 1 (p), the standard error and breadth of 95% CIs for the estimate of p for that level is high; if all values of that level are zero the SE is very high and the 95% CI ranges from 0 to 1. As a result, pairwise post-hoc comparisons made with this level are rarely significant (or never if all values are 0), even though I might expect them to be. I find this odd because I thought that standard error for a binomially distributed variable becomes smaller when, for a given n, p approaches 0 (or 1).   I have posted some example data and code for a simple example with one independent categorical variable (Cat1). For level "RE", all values but one of the responding variable "Y1" are 0. For the responding variable "Y2" all values for level "RE" are 0. If you run the code, below, you can see that the SE for RE is high in the former case, and very high (with 95%CIs of p ranging from 0 to 1) in the latter case.   Note-I am using GLIMMIX because I am ultimately included random factors in this model, but I did a bit of experimenting in PROC genmod and proc logistic and I think the same thing more or less occurs.   Can anyone explain what causes this phenomenon, or -in the case it is an error in my coding-what I may be doing wrong?   Thanks for any help!    *case where Y1 has only one non-zero value in Cat1="RE";   PROC GLIMMIX data=Fake_Data2; CLASS Cat1 ; MODEL Y1 = Cat1 / DIST=bin LINK=logit SOLUTION; LSMEANS Cat1 / CL ADJUST= tukey ilink; RUN;   *case where Y2 has all zero values in Cat1="RE";   PROC GLIMMIX data=Fake_Data2; CLASS Cat1 ; MODEL Y2 = Cat1 / DIST=bin LINK=logit SOLUTION ; LSMEANS Cat1 / CL ADJUST= tukey ilink; RUN; ... View more

Hello, I am performing analyses on water quality data. I want to test whether there are difference in average TP concentrations between 5 sites, and whether these differences depend on the year of sampling. I have three years of data: 12 samples were taken between March and October in 2014, 11 samples were taken between March and October in 2015, and 9 samples were taken between March and October in 2016. For every sampling date measurements were taken at all 5 sites. I have attached some “fake” data so you can see how it is structured. I want to “pair” samples by date to account for variability caused by weather etc. on a given date. I also want to “correct for” temporal autocorrelation- samples taken close together in time are more similar than those taken further apart. I am assuming samples taken in a given year are independent from samples taken in other years. The intervals between measurements are not equal so I am theoretically unable to use variance-covariance structures such as AR(1), toep, arma(1,1), and ante(1). I have attempted the above in proc mixed with the code at the bottom of the post. To block by day I have included a repeated statement with Date(Year) as the subject. Analyses without accounting for temporal autocorrelation showed an unstructured variance-covariance structure to be the best fit. I am allowing different var-cov matrices for each year with the group=year option. To account for temporal (residual) autocorrelation I added another repeated statement, with Site as the subject, that models the structure of temporal autocorrelation using the sp(exp) option and the variable “Day_nmbr”, which is a continuous variable representing the day as a number between 1 and 365. Note that I also have “Date” as a class variable for purpose of blocking, above. Again, I used the group=year option to allow a different autocorrelation structure in each year. Note-I haven’t explored the best residual autocorrelation model yet, but sp(exp) is a likely fit. When I run this code with either one of the repeated statements it works, but when I run it with both statements SAS only uses the second statement (I get this message: WARNING: Only the last REPEATED statement is used). I am presuming (or perhaps hoping!) that I have just coded something incorrectly. Perhaps there is a way to combine the repeated statements? Can anyone help me think of a way to code that will allow me to do what I have described? Thanks for any help!!   Code:   Proc mixed data=sites; class Date Year Site; model TP = Site Year Site*Year / solution ddfm=kenwardroger outp= pred_data; Repeated / subject= Date(Year) type=un group=year r rcorr; Repeated / subject= Site type=sp(exp)(Day_nmbr)group=year r rcorr; run;   I am using SAS 9.3 ... View more

Hello, I was wondering if anyone could clarify the use of deviance and/or scaled deviance to assess model fit (in Proc Genmod). I understand that a deviance/df value that is much greater or smaller than 1 could indicate over or underdispersion of the response variable (or model misspecification). This could lead to incorrect estimation of the standard error of parameters, and thus misinterpreation of thir statistical significance. However, it is also my understanding that with the use of certain error distributions, such as the gamma or negative binomial, a scale parameter is estimated (I am assuming to model this over/underdispersion?). In these cases, is this parameter used to correct for the effect of over/underdispersion on error estimates, i.e. does SAS' output give estimates of parameter SE that are corrected for over/underdispersion? How does this translate to affect "LR Statistics For Type 3 Analysis", if it does at all? Along similar lines, if the scale parameter is used to explicitly model dispersion, I presume scaled deviance/df can be used to assess model fit rather than deviance/df? I am asking these questions in the context of the following example, which left me unsure of the appropriate "next step": Variables: -Highly right-skewed continuous response variable that becomes normally distributed when log transformed -1 continuous independent variable, a couple categorical independent variables, select interaction terms Modelling Approaches: -general linear model approach with with log-trasformed response variable had dev/df=3.7 (which I thought was odd given the normally-distributed response variable, and the fact that the residuals were also normally distributed) -generalized linear model approach with log link function and gamma error distribution had dev/df=2.9, scaled dev/df=1.3, scale=0.45 The reported dev/df value made me question the validity of the chosen error distribution and link function, but I couldn't think of more appropriate ones given the circumstances. I am not sure where to "go" next, I hope this doesn't warrant nonparametric analysis!!!! However, if the scaled dev/df value is the appropriate term  to use to assess model fit, I think I am okay... Any clarification on any of the above questions/musings would be great!! Thanks! Madison ... View more

Hello, I am performing an analysis that is kind of new to me, and is more complicated than any examples I can find online. I am a little bit over my head, and am experiencing convergence issues when I include certain interaction terms... I was wondering if anyone could answer some questions... DESCRIPTION -Goal is to test the difference in nitrogen export ("TN_Exp") from watersheds with differect agricultural intensities ("AgInt", 3 levels) -measured monthy (March through Oct) exports of nitrogen in 19 watersheds, for 12 years -Also wanted to test the effect of monthly precipitation ("RunnoffDepth", covariate) on export -And ecological area ("EcoA", 4 levels, but not fully crossed with AgInt, so I have nested EcoA within AgInt: "EcoA(AgInt)") -And whether there was a linear trend in the 12 years ("Year", covariate) -And whether there was a seasonal effect ("Month", 8 levels) -As this is an RM design, I Included an r-side random  month term with Stream(AgInt) as a subject. I didn't necessarily expect temporal dependance between years, aside from any linear trend, however, see questions... -TN_Exp is a right skewed continuos variable, thus I used a gamma distribution, with a log link function -Stream(AgInt) was included as a random effect, as I did not care about specific differences between streams, but wanted to control for unmeasured characterisitcs at the watershed level CODE Proc glimmix data= ExpCoeff; Class AgInt Month Stream EcoA; Model TN_Exp =   Month  AgInt  EcoA(AgInt)    RunoffDepth  Year   AgInt*Month   RunoffDepth*AgInt  RunoffDepth*Month  Year*AgInt  Year*Month   RunoffDepth*Year / dist=gamma link=log solution ddfm=kenwardroger; Random  Stream(AgInt) / solution; Random  Month /  Subject= Stream(AgInt) Type= cs residual ; covtest 'indep' indep; nloptions maxiter=2000 tech= newrap; Run; QUESTIONS 1) In the "solution for random effects", all of the p-values for the estimated coefficients of each stream are >0.05 (in fact, the lowest is 0.41) Does this mean that the random effect of stream can be removed? 2) Is it necessary to nest Stream within Agint? Each stream has a different name...or perhaps I should be doing more complicated nesting since there are mulitple categorial variables???? 3) For the repeated statement, should month be nested within year...I don't think so as I am assuming the correlations between months are similar from year to year... 4) For the Random Month effect, I would expect an Ante(1) or TOEP variance-covariance matirx would be better than CS, but get the following error when I use anything more complicated than cs: WARNING: Obtaining minimum variance quadratic unbiased estimates as starting values for the covariance parameters failed. How might I be able to fix this? 5) If I wanted to consider temporal autocorrelation between years, would I include an an ADDITIONAL random statement?? i.e. Random Year / Subject= Stream(AgInt) 5) When I include interaction terms with EcoA(AgInt), the model does not converge... I am celarly not understanding something... can anyone explain this or enlighten me? Perhaps this is becasue I have interaction terms with AgInt, within which EcoA is nested? 7) When I include three-way interaction terms (i.e. RunoffDepth*AgInt*Year RunoffDepth*Month*Year Year*AgInt*Month) the model does not converge. Should I  allow more iterations, or is the model perhaps simply too complicated. Is it acceptable to not include three-way interaction terms? It would be onerous to interpret such a complicated model anyways.... Is it acceptable to simply choose interaction terms you have reasons to believe would exist/ you are interested in? 😎 If you happen to notice an potential errors, I would love to know!!! THANK YOU SO MUCH TO ANYONE WHO HAS ANY INSIGHT!!!! Madison ... View more