Your model in your second post specifies a 3-way factorial in a randomized block design, but in your first post you described your design as a split-plot with (I think, I'm not entirely sure) whole plots in blocks. So this second model specification does not appear to match the design. Your model in your first post corresponds (I think) to your actual design.
You can do a lot of exploration and testing of lsmeans within the LSMEANS statement in either the MIXED or GLIMMIX procedure. I've illustrated some of these options in my code below; refer to the documentation for more details about what these options do. Using the DIFF option on the LSMEANS statement does not get you very far in interpretation of interactions; as I noted before, a thoughtful, context-specific approach is better. I've used GLIMMIX rather than MIXED to get interaction plots and to use SLICEDIFF, but otherwise the LSMEANS options I've used are available in MIXED.
data test;
input Rep Tillage$ Crop_Residue$ N_Dose$ SOC;
datalines;
1 Contour with N0PK 30256.7
2 Contour with N0PK 16870.24
3 Contour with N0PK 26706.44
4 Contour with N0PK 23508.73
1 Contour with N120PK 24466.53
2 Contour with N120PK 24032.51
3 Contour with N120PK 26145.86
4 Contour with N120PK 23643.35
1 Contour with N60PK 27335.54
2 Contour with N60PK 23438.42
3 Contour with N60PK 27678.23
4 Contour with N60PK 26230.67
1 Contour without N0PK 23386.71
2 Contour without N0PK 18918.63
3 Contour without N0PK 27250.83
4 Contour without N0PK 23500.25
1 Contour without N120PK 24153.19
2 Contour without N120PK 23274.78
3 Contour without N120PK 23882.67
4 Contour without N120PK 21778.33
1 Contour without N60PK 27734.6
2 Contour without N60PK 19916.57
3 Contour without N60PK 25009.3
4 Contour without N60PK 22351.28
1 Reduced with N0PK 23766.65
2 Reduced with N0PK 30368.58
3 Reduced with N0PK 29149.02
4 Reduced with N0PK 27259.67
1 Reduced with N120PK 28545.35
2 Reduced with N120PK 28205.5
3 Reduced with N120PK 29953.75
4 Reduced with N120PK 28287.69
1 Reduced with N60PK 32671.92
2 Reduced with N60PK 25906.4
3 Reduced with N60PK 25047.6
4 Reduced with N60PK 27056.77
1 Reduced without N0PK 25677.25
2 Reduced without N0PK 23565.46
3 Reduced without N0PK 26469.71
4 Reduced without N0PK 28219.43
1 Reduced without N120PK 29238.49
2 Reduced without N120PK 24391.22
3 Reduced without N120PK 28182.27
4 Reduced without N120PK 27081.03
1 Reduced without N60PK 28364.79
2 Reduced without N60PK 27046.47
3 Reduced without N60PK 25500.1
4 Reduced without N60PK 24558.67
run;
/* Rename n_dose levels to be quantitative */
data test;
set test;
if n_dose = "N0PK" then n_level = 0;
else if n_dose = "N120PK" then n_level = 120;
else if n_dose = "N60PK" then n_level = 60;
run;
/* Plot observed data, for visual assessment of various things */
proc sort data=test out=test_out;
by n_level;
proc sgpanel data=test_out;
panelby tillage crop_residue / columns=2;
series x=n_level y=soc / group=rep markers;
run;
/* ASSUMING a 3-way full factorial treatment structure
in a split-plot design with whole plots in blocks */
proc glimmix data=test
plots=(studentpanel boxplot(fixed));
class rep tillage crop_residue n_dose;
model soc = tillage | crop_residue | n_dose;
random rep /* block variance */
rep*tillage ; /* whole plot variance */
/* residual is subplot variance */
lsmeans tillage | crop_residue | n_dose;
/* PRETEND that tillage*n_dose interaction is significant */
lsmeans tillage*n_dose / plot=meanplot(sliceby=tillage join cl) /* interaction plot, not available in Proc MIXED */
slice=(tillage n_dose) /* simple effects tests */
slicediff=tillage adjust=tukey; /* pairwise comparisons among n_dose means adjusted for Type I error within each level of tillage */
run;
I hope this is of help.
... View more