06-01-2015 01:59 PM
I am doing a smiulation study on analysis of spatially correlated data. My response conditioned on random components has a Poisson distribution with different means for 5 different treatments. Also the response has a spherical spatial structure. Since the response has different means depending on treatments I use GROUP=trt option in the RANDOM statement. Then the expected variance of the response will have 5 different values each for each treatment because variance depends on mean. However, the expected range of the response is same which is 6 for all the treatments. So when I use PARMS statement to specify initials for covariance parameters, I should give 5 values for variance components and one value for the ranges of all the treatments. And I don't know how to do this in GLIMMIX. There might be other ways such as writing two RANDOM statements one for variance and one for range. Again, I don't know if that works with the TYPE=sp(sph).
Here is more detail about my work. The Poisson means for 5 different treatments are 10, 10.8, 11.6, 12.4, and 13.2. So the initial values of variance components are reciprocals of the Poisson means. There are 5 plots received each treatment and 9 sub sampling points within each plot. To specify spatial structure among subsampling points, the variable pos is used. Y is my response. So the GLIMMIX code that I have been using is
CLASS trt plot pos;
MODEL Y=trt/ DIST=Poisson DDFM=KR;
RANDOM pos/ SUBJECT=intercept GROUP=trt TYPE= SP(SPH)(row col);
PARMS (0.1 6 0.09259 6 0.08621 6 0.08064 6 0.07576 6);
Notice that in the above PARMS statement, I have given reciprocal of Poisson mean followed by the range 6 for each treatment. When I use the above code I get different estiamtes for the range of response for different treatments. However, I want to get a single estiamte for the range of the response.
I will appreciate if anyone could help me to figure this out in GLIMMIX.
Thank you in advance.
06-02-2015 07:23 AM
Before using the PARMS statement, what sort of results do you get when the procedure uses the default starting values? From that, you may get better estimates of what goes into the PARMS statement.
Also, the subject= should be the plot (i.e., subject=plot) rather than 'intercept' which is undefined as a variable. I am not sure what the variable 'pos' indicates--I am inferring it is the position of the plot, which is dependent on the row and col variables. I would rewrite the random statement as:
RANDOM intercept/ SUBJECT=plot TYPE= SP(SPH)(row col);
removing the GROUP= option. The use of DIST=Poisson forces the means and variances to be equal, and it may not be necessary to add this in. I would run both ways and compare results.
06-02-2015 11:22 PM
Thank you for your reply, Steve.
When I run without PARMS statement, I get pretty large values for range and values less than 2 for variance. Since I have simulated these data, I expect the range to be 6. Also, based on GLMM theory, the variance of conditional Poisson variable is the reciprocal of the Poisson mean. Since the response is a combination of Poisson counts with different means, the estimated variance varies a lot.
Here is more information about the experimental design to under stand what is "pos"and what goes to SUBJECT. The field consists of 15 by 15 grid so that there are 225 points. Each plot is 3 by 3 and there are 25 plots. 5 treatments are randomly assigned to plots so that 5 plots receive a single treatment. Then there are 9 subsampling points in each plot. The variables row and col specify row and column, respectively, of subsampling points. The variable pos specifies the subsabmpling point which goes from 1 to 225. So I would still think SUBJECT=pos should be correct.
I included GROUP=trt option there beacuse variance of response is different by treatment. If that option is not used, then the variance of the response deviates from the expected value due to different Poisson means within the response.
Since I have generated data using Poisson distribution, I should use the DIST=Poisson option. And the goal of this simulation study is to analyze Poisson counts.
Any more help is appreciated.
06-03-2015 09:44 AM
The variance of the conditional Poisson is the mean, not the inverse of the mean (using standard parameterization). You can use sub=intercept when you want all observations to be correlated, not just within plots (special syntax to show that there is one big subject).
You are modeling spatial covariance as a G-side effect, not an R-side one. Thus, the variance of the G side term is not related to the conditional variance of the Poisson. These are different things.