Regarding "Indeed, time spent is just one of the response we measured. there is another response value for the preference, i.e., LEDnum (num of birds under LED at night, as birds won't move at night, so we count number of birds in each cage, we have total 3 birds for each group). for time spent or proportion of time, the value is continous data, how about the number of birds (0,1,2, or 3) or the corresponding proportion? Can we use the same model to analyze it?":   The analysis of LEDnum should use a binomial distribution, which models proportion as number of "successes" out of number of "trials". "Success" is a bird under LED. Number of trials is 3. The beauty of the GLIMMIX procedure is that it offers different distributions (like beta and binomial), so to some extent the code is similar. Unfortunately, there are lots of ways that a binomial mixed model can go wrong: estimation problems, overdispersion, etc.    I don't have time to do much more right now, but see https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_a0000001403.htm https://support.sas.com/resources/papers/proceedings15/SAS1919-2015.pdf http://www.lexjansen.com/wuss/2006/Analytics/ANL-Dai.pdf   ... View more

Ah, good.   Another thought about the day effect: you are hoping/assuming that birds do not develop an aversion to one cage relative to the other (which is not the same as preferring one light type over another). If you were, for example, shocking birds when they were in the right-side cage, they might sensibly decide not to go there in the future, regardless of whether it had their favorite light bulb. This would be an example of carry-over effect that can be a problem in Latin square designs.   Regarding BirdType: Yes, it can be added to the model. Pair and BirdType are treatment factors that are associated with a Group of birds (i.e., a Group is the experimental unit for Pair and BirdType). If each trial is 10 h long, then you may be running only 2 Groups each day, requiring 24 calendar-days to complete trials for 12 Groups run 4 times. (Using "calendar-day" I mean to distinguish between the time that you ran a trial (e.g., February 15) from the "day" that marks the observation within a trial (i.e., day 1, 2, 3, 4.) For 3 different types of Group, the experiment requires 72 calendar-days. That is a fair amount of time passing; depending on the environmental conditions, bird "history", how different BirdTypes were allocated to calendar-days, and such, there could be confounding factors associated with time passing. Give that some thought, discuss the possibility with your colleagues.   proc glimmix data=have method=laplace; class pair group day birdtype; model proportionLED = pair | birdtype | day / dist=beta intercept s; random group*pair*birdtype; estimate "Intercept" intercept 1 / ilink cl; lsmeans birdtype / ilink cl; run; ... View more

Let's meet in the middle with the RANDOM statement. Try your restructured dataset with this model  proc glimmix data=have method=laplace; class pair group day; model proportionLED = pair day pair*day / dist=beta intercept s; random group*pair; estimate "Intercept" intercept 1 / ilink cl; run; This RANDOM statement estimates a variance among the 12 groups of birds.    We hope that there is no statistical evidence of a pair effect, in other words, we hope that the two pairs are essentially identical. We really hope that there is no pair*day interaction; a significant interaction might imply a difference between the two cages within a pair and/or an effect of the cage assignment to light bulb type. If there is a day effect, then the pattern of the day means might imply something like "learning" by birds over time.   If all of the fixed effects are not significant, should you omit them from the model? That's arguable. Certainly they are "nuisance" factors: they are not components of your research hypotheses. If not significant, then retaining or omitting likely will have little effect on the Intercept estimate. But you never know.   This model does not address whether the covariances among the four repeated measures are unequal (it's essentially invoking compound symmetry). The beta distribution does not lend itself to explicitly modeling of the covariance structure of repeated measures; as I understand it, there is a conflict between the estimation of the scale parameter for the beta distribution and the estimation of covariance parameters for repeated measures.   This model does not address the assignment of light bulb type to cage within a pair that changes over time--the design is a Latin square, but not a simple one. These design elements might be added, but they would complicate the model appreciably. If there is no statistical support for pair*day, I probably would not try to model the Latin square aspects.   Because I've added fixed effects to the model, I would expect 10 denom df for intercept.   Let me know what you get!   Was each trial for a group of birds the same length of time? In other words, if you add time spent under LED and time spent under CFL, do you get the same value for each trial?   Although the model provides an Intercept estimate (and CI, etc.) for just one light type, you can get estimates for the other, either algebraically or by analyzing proportionCFL instead of proportionLED. ... View more

  I do not understand the structure of your dataset and how it reflects the design of your experiment you described in your initial question. What is lightdur? What is lightpre? You say you have 12 groups of birds; how are these 12 identified in your dataset?   Did you repeat the trial 4 times (for example, once each day for four days) on the each group of birds?   I suspect the values of lightdur and lightpre are means over groups of birds. If so, then you cannot do statistical analysis on this dataset. For statistical analysis, you need data for independent replications (which I think would be groups of birds).   One format for a dataset that I would expect has 48 observations, one for each day for each group of birds. Variables would be (1) group identification (with 12 unique values); (2) pair (A or B); (3) day (1, 2, 3, or 4); (4) amount of time spent under LED; and (5) amount of time spent under CFL. Subsequently in a datastep, you could compute the proportion of the total time (amount under LED + amount under CFL) that a group spent under LED; this proportion is the response variable that you would use in the model I proposed.   Cage (the two areas within each pair) could be incorporated into the model, but it would make for a more complex statistical model, and we are already challenged by this simpler model. I would not include it at this stage. I might consider including it later if there appears to be a day effect of some sort.   If you want to assess preference for one light or the other, then the proportion of time spent under one light is a good metric. If you agree with that, then think about how the CFL proportion is just another way to measure the LED proportion:    CFL proportion = 1 - LED proportion   Analyze one or the other, but you do not need both because the two metrics are just different ways to measure exactly the same thing.   If you want to assess the amount of time spent under one light or the other, well, that could be another analysis and is not the objective you described in your initial question.  ... View more

You're welcome.   First comment: Note that your RANDOM statement is not what I recommended. The different RANDOM statements generate different models with different assumptions about the role of variables in the design.   Second comment: The statement estimate "Intercept" intercept 1 / ilink cl; will provide an "Estimate" and "Standard Error" which are on the link scale; here, the link is logit. The columns labeled "Mean" and "Standard Error Mean" (generated by the ilink option; ilink means inverse link) are estimates on the proportion scale (i.e., the proportion of time spent under LED light), something (but not exactly) like back-transformations of "Estimate" and "Standard Error", respectively. "Estimate" will not look anything like a mean proportion computed from the raw data. "Mean" will be close, but it will not be the same as the mean computed from the raw data for various reasons that are too complicated to describe here--but are covered in Stroup (2013) Generalized Linear Mixed Models, CRC Press.   The paired t-test is wrong. The proportion of time under LED is not just (negatively) correlated with the proportion of time under CFL; these two metrics are totally redundant--there is no information in one metric that is not already in the other. If you told me that in one trial, the proportion of time under LED was 0.62 (62%), I could tell you that the proportion of time under CFL in that trial was 0.38, because these two values must sum to 1 (100%). Applying a t-test to these "paired" data is not valid. Don't do it. ... View more

On the number of parameters topic, I disagree. I'm happier when I have a data set to test, so I dummied one up (very badly; if I had time, I'd pull down my copy of Rick Wicklin's Simulating Data with SAS and figure out how to simulate data with specified covariance parameters, to check against estimates). [Edit: For the model I posed earlier, the number of random effects parameters is 12; see the attached code example.] I'm thinking my suggested model (or a variant) would work--if it would converge with kpboh's data, which is clearly an unresolved problem. See what you think about the code in the attached file.   ... View more

(2) The documentation for MIXED | Details | Computational Issues | Convergence Problems includes this statement:   "If PROC MIXED stops because of an infinite likelihood, recheck your model to make sure that no observations from the same subject are producing identical rows in or and that you have enough data to estimate the particular covariance structure you have selected."   Your TABULATE results show that the problem is not due to coding (observations from the same subject producing identical rows in or ), which means that it is probably due to not enough data. Adjusting parameters that control the estimation process (covered in the documentation and the Kiernan et al. paper that @SteveDenham recommended) might help, or not.    So....what next? Other than give up. I'd follow Steve's lead and try simpler covariance structures, building up to more complex, and see how close you can get to what you want.   A split-plot-in-time-like model that specifies 3 sources of variance (in addition to residual variance) with default covariance structures /*RandomStmt1*/ random intercept / subject=pop; /*RandomStmt2*/ random intercept / subject=id(pop); /*RandomStmt3*/ random var / subject=pop; Try 1: RandomStmt1 and RandomStmt2 (i.e., drop RandomStmt3) Try 2: RandomStmt1 and RandomStmt2 and RandomStmt3 Try 3: Steve's suggestion: RandomStmt1 with group=var Try 4: RandomStmt1 with type=un or unr or fa0(3) (see http://comp.soft-sys.sas.narkive.com/RqRMBFgu/difference-between-un-vs-unr-type-option-in-random-statement) Try 5: RandomStmt2 with default covariance structure, and RandomStmt3 with type=un or unr or fa0(3) Try other combinations that you think are appropriate, depending on whether simpler Trys converge.   Clearly, I haven't attempted these, and so if one looks wrong, it might be wrong.   None of these Trys are doing anything fancy with residual variance. Should they? You don't appear to be interested in covariance structure for residual variance--your focus is on covariance structure for pop variance--and maybe you can ignore residual variance to some extent, although I'm not sure that's statistically valid.     ... View more

Hi Steve @SteveDenham and @kpboh,   Thoughts:   (1) Regarding the number of parameters: 78 seems too high. I was thinking it would be 6 for the first random, 6 for the second random, and another 3 for the var means. But 15 is still a lot, and maybe I'm not counting correctly. I thought about taking this problem for a spin with simulated data, but I haven't had the time.   (2) kpboh, as Steve suggests, you can check the data structure. I use the TABULATE procedure a lot for this task:   proc tabulate data=tall; class pop id var; table pop*id, var; run; I think you would hope to see "1" in each cell.      (3) kpboh, in your original message, your example dataset was "pops", but the dataset used in GLIMMIX was "tall". Are these two datasets equivalent?   (4) kpboh, it's not clear whether you ran a model with both class var pop id; and repeated var /subject=id(pop) type=unr; and still had lack of convergence. You might have, it's just not clear.   (5) If my model suggestion is not feasible and if you want covariances at the pop level, could you average the y1, y2, and y3 values over the ids within each pop? Do you have the same number of ids in each pop? (Does it even make sense to average over ids, in context or statistically?) How else can we get at covariances at the pop level....? (There's the rub.) proc mixed data=tall_MEANS method=reml covtest asycov; class var pop; model y=var ; repeated var / subject=pop type=unr; run;   Susan   ... View more

@SteveDenham, following your lead (which looks like a good direction to me), what do you think about something like this: proc mixed data=tall method=reml covtest asycov; class var pop id; model y=var ; random var / subject=pop type=unr; repeated var / subject=id(pop) type=unr; run; It's a bit manic.  Does the covariance structure of the entire model make sense? I'm not sure.   I added id to the class statement which I think will help with the estimation problem. I dropped noint from the model statement but doubt that matters. Would there be any hope of convergence with only 13 levels of pop (edit: and possibly few levels of id within pop) and so many covariance parameter estimates?   ... View more

Thank you for the informative responses. My turn!   1. Because of the censoring and because the body condition covariates are time-varying, I'm not seeing how a good binary GLMM can be constructed. The resume values don't switch back and forth across periods between N and Y depending on body conditions; a resume value is N until it (possibly) becomes Y. Looks like event data to me.   2 and 3. Your approach to variable selection is arguably treacherous. But I concur that you have too many predictor variables for your sample size. Even with this process of down-sizing to 7 variables, you still need to assess multicollinearity among them. A quick Google of "multicollinearity detection" generates a lot of useful links. Multicollinearity is a property of the predictor variables, so you can read about multicollinearity in regression texts that deal with normal-distribution response and use multicollinearity statistics generated by regression software (like REG).   4. If body conditions are largely a consequence of intervention, then it would be appropriate to drop intervention.    5. By default, if an observation contains a missing value for one or more predictor variables, the software will drop that observation in the model fitting; this is true for most (if not all) software packages. So, you need to deal with those missing values in some fashion, or you'll lose those observations.   6 and 7. I commend you for not blinding using stat software! This forum is usually good to help with code when you provide enough detail to identify an appropriate model for a corresponding dataset. In this case, as you can tell, I don't believe a binary GLMM is appropriate so I cannot help with syntax for that model (or should not, but see below). If you switch to a time-to-event model, this website which I linked to earlier http://www.ats.ucla.edu/stat/r/examples/alda/ has SAS code examples, and the forum could help with code if you hit snags.   In the interest of your ongoing education, I will add that if I liked your approach, which I don't, I would start with     proc glimmix data=resumption method=laplace; class period id; model resume = wt bmi perbf fm fmi lmi lpct tpct tlpctr apct gpct agpctr period / link=logit dist=binomial or solution; random int / subject=id; random period / subject=id type=ar(1); run; This sort of model almost always requires some twiddling, even beyond that needed to identify a good covariance structure type. Including "type=ar(1)" and "residual" on the first random statement is wrong. Including both random statements as here identifies an AR(1)+RE covariance structure. Omitting the first random statement and retaining only the second identifies an AR(1) covariance structure. See Littell et al. http://onlinelibrary.wiley.com/doi/10.1002/1097-0258(20000715)19:13%3C1793::AID-SIM482%3E3.0.CO;2-Q/abstract  or Stroup (2013) Generalized Linear Mixed Models for details about the distinction between AR(1)+RE and AR(1). If you had enough subjects and enough periods and more complete longitudinal data, you could consider a model with random slopes, specified as       random wt bmi perbf fm fmi lmi lpct tpct tlpctr apct gpct agpctr / subject=id type=vc; but this is way too complex for what you have to work with.     ... View more

I'll offer some thoughts. Hopefully other people will weigh in as well.     1. Consider a "time to event" analysis, where event is resumption of menses. Given the extent of censored observations (i.e., women who withdrew, and women who did not resume menses before the end of the study), I doubt that a GLMM with a binary response will work well. The text by Singer and Willett (Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence) and the code here http://www.ats.ucla.edu/stat/r/examples/alda/ could be a good place to start.   2. 30 subjects, about half of whom exhibited the event, is a small sample. Keep in mind your effective sample size and the risk of over-fitting.   3. With correlations that high among the body condition covariates, you will surely have issues with (multi)collinearity. Consider a more complete assessment of multicollinearity, and then some form of dimension reduction, e.g., dropping predictor variables, forming an index, principal components analysis, factor analysis. Again, keep in mind your effective sample size and the risk of over-fitting.   4. You say you are not interested in the "intervention" effect. Are different women exposed to different interventions in your dataset? If so, I would not think that you could ignore it unless intervention absolutely does not matter.   5. The missing covariate values will pose problems. Consider imputation? Or drop those covariates with missing values if they are highly correlated (and thus redundant) to other covariates.   6. These thoughts are all independent of the software that you use, whether SPSS or SAS (or anything else).   7. I suspect this analysis will be a challenge, quite probably more than can be dealt with adequately in this forum. Consider finding someone at your institution with statistical expertise to guide you.   ... View more

The dataset included in your code has 8 subjects. Is this the full dataset or just a subset? If it's a subset, how many subjects do you have?   ... View more