Thanks for pointing this out. I did not realize what I did wrong. Then, if I change the code to this form, will I be able to use the residual output from proc glimmix?   Thanks in advance.   proc glimmix data=a1.prediss1 plots=(studentpanel boxplot(fixed student)); class phase_info1 Phoneme_stress subject; model GMP_log= phase_info1/solution dist=gamma; output out=R pred=pred resid=r; random intercept phase_info1 / subject=subject; random _residual_ / group=phase_info1; run; ... View more

Yes. I generated residual variable after proc mixed and changed the model to proc glimmix with gamma distribution. I expected the resultant residual should be different, but the normality test results were the same. I questioned whether the code was correct. How about the codes afterwards to test the normality assumption, do you see any errors? Thanks in advance.   proc glimmix data=a1.prediss1 plots=(studentpanel boxplot(fixed student)); class phase_info1 Phoneme_stress subject; model GMP_log= phase_info1/solution dist=gamma; output out=a1.prediss1 pred=pred resid=r; random intercept phase_info1 / subject=subject; random _residual_ / group=phase_info1; run; data R1; set R; absr=ABS(RESID); run; proc univariate data=R1 normal; BY phase_info1; var Resid; histogram / normal ctext = blue; run; proc sort data=R1; by phase_info1 Phoneme_stress subject; run; ods graphics on; proc univariate data=R1 normal; var GMP_log; histogram / normal ctext = blue; run; ... View more

I can only attach one article at a time. Here is the second article. ... View more

Hello. I have been helped by Sid before. Hello, Sid.   I may not be a good position to advise anyone in this community, but I can share good sources you can utilize to have a better understanding about mixed model and how you link the model to the SAS codes. If you read the attached two articles, it will help you a lot when designing your mixed model. I think the most important thing is to determine what will be your level-1 variable and what will be your level-2 variable, etc.  If you have more than two levels, it gets very complicated though. Once you are sure about your equation, you know what goes to fixed and random statements.   I think that is why Sid asks what is the relationship among year, rep, and trt: whether the trt 1 on the year 2015 is the same trt 1 in year 2016, etc.   Hope this helps.             ... View more

Thanks for all of your help and advice. Your reply is like quenching my thirst.   Here is what I developed so far with a help of others. To run the proc glimmix and to get the R dataset for homoscedasticity testing, I modified your code to this:   proc glimmix data=dissertation1 plots=(studentpanel boxplot(fixed student)); class group slc subjectid trial1; model ISI_Avg_sqrt = group|slc / ddfm=kr2; output out=preddata pred=pred resid=r; random intercept / subject=subjectid(group); random slc / subject=subjectid(group); random intercept / subject=trial1(subjectid group); /* because I have not yet confirmed about homogeneity of variance*/ /*random _residual_ / group=slc;*/ lsmeans group*slc / plot=meanplot(sliceby=group join cl); slice group*slc / sliceby=group diff alpha=.0083; /*lsmestimate group*slc "Group effect: SLC 1 v 2" 1 -1 0 -1 1 0, "Group effect: SLC 1 v 3" 1 0 -1 -1 0 1, "Group effect: SLC 2 v 3" 0 1 -1 0-1 1 / adjust=simulate(seed=29847);*/ run;     Or to use this: proc mixed data=dissertation1 ic; class group slc subjectid trial1; model ISI_Avg_sqrt= group|slc /ddfm=kr2 outp=R; random intercept / subject=subjectid(group); random slc / subject=subjectid(group); random intercept / subject=trial1(subjectid group); repeated / group=slc; run;   /*Then, test homoscedasticity assumption*/ data R1; set R; absr=ABS(RESID); run; proc glm data=R1; class group slc subjectid trial1; model ABSR=group|slc; run;   However, I still have a remaining question regarding how I will deal with "trial" variable. Although I did not include any research question asking the effect of trial, I think when I planned my statistical analysis models, I included "trial" as part of my model. So my question is whether above proc glimmix or proc mixed codes reflected what I planned. I cannot figure out what goes where in SAS codes to reflect above mixed model equation more authentically. There is also possibility that I can change my equations if needed.   Thanks in advance for your help.       ... View more

If I readdress my question, could you please help me understand following two lines?   random intercept / subject=trial1(subjectid group);   (--> Does it mean to report the variance due to a variable "trial"? Is there a possibility that "trial" was treated as "repeated" variable in proc glimmix through this code, because proc glimmix does not have repeated statement?) random _residual_ / group=slc;                         (--> I think I have a vague sense that residuals are divided by the level of slc, but I am not quite understanding why we divide them. How does it take care of heteroscedasticity problem?)   Thanks in advance. ... View more

Sorry for my late reply to this comment. I realized you suggested pertinent contrast, however it is not the term I know.   Your codes looks like it provides the marginal effect, but not the simple main effect at each line of the other dependent variable.   I have a flow chart to follow when interaction is significant or not significant. In that flow chart, it was recommended to examine interaction contrast or simple main effect if the interaction is significant. That was why I looked for the codes how I can examine simple main effects. If you can teach me what the pertinent contrast is and why it is better than pairwise comparison, I will use the information your codes provided, but for now because I am not fully understanding what you are suggesting, I think I will go with the direction what I have learned at classes.   Thanks for your suggestion though.       ... View more

I think the dataset I uploaded in my previous reply contained some errors. I uploaded the right one again.   Assuming the current code I described in my previous reply is correct, I start to write up the results.   Now I wonder how I write the codes to examine posthoc analysis when interaction effect is significant. For example, if I want to come up with p-value for the effect of one variable (e.g., Trial) at each level of another variable (e.g., SLC). Will these two lines for such simple main effects?   lsmeans trial1 / diff adjust=simulate(seed=98375) lines; lsmeans group*slc / plot=meanplot(sliceby=group join cl);    Thanks for your help again in advance.   ... View more

Why did you have each participant do 10 trials? What research questions do you want to answer about trials? A: It is to obtain enough number of trials for each target sequence. Just one trial may not work.   Is SLC the sequence length condition (3-, 6-, or 9-syllable)? A: Yes.   Did you use the same order of sequences in all 10 trials? A: I am not sure what you are asking here. The variable "Trial1" reflects the order of 10 trials that were produced for each sequence . So the trial 1 was produced before trial2, trial 2 was produced before trial 3, etc. This order does not necessarily match with the numbers in the file name though. The 12 sequences were randomly mixed up in each block. Block 1 had different order of 12 sequences from block 2, the block 2 had different order of 12 sequences from block 3, etc.   And, to confirm, the same order of sequences with all participants? In other words, there was only one order of sequences for the whole study? A: Yes.   Do you think that order of sequences would affect the response? A: No. That is why I randomly mixed up but made sure 12 sequences were produced twice in one block.   How did you choose which 3 sequences of the 12 to use for analysis? A: Randomly.   I think the dataset I uploaded in my previous reply contained some errors. I uploaded the right one again. ... View more

The participants produced 12 sequences (4 3-syllable sequences, 4 6-syllable sequences, 4 9-syllable sequences) 2 times in one block. They produced these five times over five blocks. They produced each sequence 10 times. So basically what I am interested in is the comparison among three sequence length conditions and blocks mean nothing. The tasks were blocked only to give participants a few minute breaks. Although my research questions do not question about trial effect, I am still interested in it and I assume there might be learning effect when participants repeat the same task over ten times (trial1). So I wanted to include this as one of random or repeated variable in the model.   Among 12 sequences, I picked 1 3-syllable sequence, 1 6-syllable sequence, and 1 9-syllable sequence as my target sequences for analysis. The other ones were fillers and were not included in the dataset I provided.   12 sequences were randomly mixed, but I pre-determined this random order before the experiment, so every participant underwent the same order of tasks.   I don't know if this design changes any of your codes. Because I face some problems whenever I include random and repeat statements together in one model and somehow your codes work fine, I am careful to change any parts of your code. So if you think your code is not reflecting my design, I am not capable of changing it without causing errors.   These are the final codes I used. I am also attaching my final dataset. Could you please help confirming this looks fine? After testing with these codes, normality assumption did not meet even with the log tranformed gmp_syl, so I decided to report only glimmix results.  Although both proc glimmix and proc mixed produced the same results.   Thanks.     /* Import data from xlsx */ libname a1 'C:\Data Analysis\Experimental Results\SAS'; proc import datafile="C:\Data Analysis\Experimental Results\SAS\dissertation.xlsx" out=Work.dissertation dbms=xlsx replace; run; /* Delete unfilled columns and unfilled rows */ data dissertation1; set dissertation; /*drop var20-var56;*/ where whole_1st= 1; run; /* Variable transformation */ data dissertation1; set dissertation1; gmp_syl_log = log(gmp_syl); gmp_syl_sqrt = sqrt(gmp_syl); run; title1 "Coding check"; proc freq data=dissertation1; table group subjectid slc block trial1; run; title1 "Pattern of observations"; proc tabulate data=dissertation1(where=(gmp_syl ne .)); class group subjectid slc trial1 block; table group*subjectid, slc*trial1 / misstext="X"; run; proc sort data=dissertation1; by group subjectid slc trial1; run; title1 "Observed data"; proc sgpanel data=dissertation1 noautolegend; panelby slc group / columns=2; series x=trial1 y=gmp_syl / group=subjectid markers; run; title1 "Observed data, log scale"; proc sgpanel data=dissertation1 noautolegend; panelby slc group / columns=2; series x=trial1 y=gmp_syl_log / group=subjectid markers; run; proc sgpanel data=dissertation1; panelby group subjectid / columns=6; series x=trial1 y=gmp_syl_log / group=slc markers; run; /*unconditional model ICC check*/ proc mixed data=dissertation1 ic; class group subjectid slc trial1; model GMP_syl_log= trial1/solution; random intercept trial1/sub=subjectid(group); run; /*homoscedasticity assumption testing*/ proc mixed data=dissertation1 ic; class group subjectid slc trial1; model GMP_syl_log= group|slc|trial1 /ddfm=kr2 outp=R; random intercept / subject=subjectid(group); repeated trial1 / subject=slc*subjectid(group) type=ar(1); run; data R1; set R; absr=ABS(RESID); run; proc glm data=R1; class group subjectid slc trial1; model ABSR=group|slc|trial1; run; /*normality assumption testing*/ proc mixed data=dissertation1 covtest method=ml;*plots=all; class group slc subjectid trial1; model gmp_syl_log = group|slc|trial1/ ddfm=kr2 outp=R; random intercept / subject=subjectid(group); repeated trial1 / subject=slc*subjectid(group) type=ar(1); /*lsmeans trial1 / diff adjust=simulate(seed=98375); lsmeans group*slc ;*/ run; proc sort data=R; by group slc subjectid trial1; run; ods graphics on; proc univariate data=R normal; BY group slc ; var Resid; histogram / normal ctext = blue; run; title1 "AR(1) among TRIAL1 levels using MIXED"; proc mixed data=dissertation1 covtest ;*plots=all; class group subjectid slc trial1; model gmp_syl_log = group|slc|trial1 / ddfm=kr2; random intercept / subject=subjectid(group); repeated trial1 / subject=slc*subjectid(group) type=ar(1); lsmeans trial1 / diff adjust=simulate(seed=98375); lsmeans group*slc ; run; title1 "AR(1) among TRIAL1 levels using GLIMMIX"; proc glimmix data=dissertation1 plots=(studentpanel boxplot(fixed student)); class group subjectid slc trial1; model gmp_syl_log = group|slc|trial1 / ddfm=kr2; random intercept / subject=subjectid(group); random trial1 / subject=slc*subjectid(group) type=ar(1) residual; lsmeans trial1 / diff adjust=simulate(seed=98375) lines; lsmeans group*slc / plot=meanplot(sliceby=group join cl); lsmestimate group*slc "Group effect: SLC 1 v 2" 1 -1 0 -1 1 0, "Group effect: SLC 1 v 3" 1 0 -1 -1 0 1, "Group effect: SLC 2 v 3" 0 1 -1 0-1 1 / adjust=simulate(seed=29847); run;     ... View more