Random graphs (48): Interaction plot with overlaid data points

// Simulate data
clear
set seed 1
set obs 500

generate  e  = 0 + (500 - 0) * runiform()     // To generate random variates over the
generate  x1 = 0 + (800 - 0) * runiform()     // interval [a,b), a+(b-a)*runiform()

generate  x2 = round(0 + (1 - 0) * runiform()) // Binary variable
generate  y  = (x1 + x2 + (x1*x2) + e) / 1000

// Calculations for Aiken & West (1991) style plot
regress y c.x1##x2

qui sum x1
local x1_minsd = r(mean) - (2*r(sd))
local x1_mean  = r(mean)
local x1_plusd = r(mean) + (2*r(sd))


// Calculate predicted values for plot
margins, at(c.x1 = (`x1_minsd' `x1_mean' `x1_plusd') ///
            c.x2 = (0 1)) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
   title("Interaction plot with overlaid data points") ///
   ytitle("y") ///  
   ylabel(, format(%6.1f)) ///
   xtitle("") ///
   plotopts(msymbol(none)) ///        // Turn off markers
   plot1opts(lpattern(longdash)) ///  // Define line types here
   plot2opts(lpattern(solid)) ///
   addplot(scatter y x1 ///
         , symbol(o) ///
           xlabel(`x1_minsd' "-2 SD" ///       // The addplot seems to override
                  `x1_mean' "Average x1" ///   // the regular axis label 
                  `x1_plusd' "+2 SD") ///      // command
           legend(subtitle(x2) ///             // And the legend command
                  order(3 "x2 = 0" 4 "x2 = 1" 2 "95 % CI"))) ///    
   legend(pos(5) ring(0)) ///                  // But not completely
   name(plot, replace)

Random graphs (27): Between and within regression

// Case A
clear
set seed 1
set obs 300

generate  e = 0 + (200 - 0) * runiform()  // To generate random variates over the
generate  x = 0 + (600 - 0) * runiform()  // interval [a,b), a+(b-a)*runiform()

generate  y = e // Relationship between X and Y random
generate id = 1
replace   y = (e + 200) if x > 200 & x <= 400
replace  id = 2         if x > 200 & x <= 400
replace   y = (e + 400) if x > 400
replace  id = 3         if x > 400

bys id: egen ym = mean(y)  // Generate cluster means for Y and X
bys id: egen xm = mean(x)

twoway (scatter y x if id == 1, msymbol(oh)) ///
       (scatter y x if id == 2, msymbol(dh)) ///
       (scatter y x if id == 3, msymbol(sh)) ///
       (scatter ym xm,          msymbol(S)) ///
       (lfit ym xm,             lpattern(solid)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
      , /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/  legend(off) ///
        xtitle("Component 1") ytitle("Component 2") ///
        title("Case A:" "Regional correlation," "individual orthogonality") ///
        xlabel(none) ylabel(none) name(casea, replace)

// Case B
clear
set seed 1

set obs 4 // Generate four clusters

generate id = _n
generate u_i = 300 if id == 1  // Generate cluster-specific intercepts
replace u_i = 0 if id == 2
replace u_i = 0 if id == 3
replace u_i = -300 if id == 4

expand 100  // Generate units in clusters

bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70)                // Generate level-1 error term

generate y = x + u_i + e_ij
drop if y < 0 | x < 0
drop if id == 2 & x >= 300
drop if id == 3 & x  < 300

bys id: egen ym = mean(y)
bys id: egen xm = mean(x)

twoway (scatter ym xm, msymbol(S)) ///
       (scatter y  x  if id == 1, msymbol(oh)) ///
       (scatter y  x  if id == 2, msymbol(dh)) ///
       (scatter y  x  if id == 3, msymbol(sh)) ///
       (scatter y  x  if id == 4, msymbol(th)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
       (lfit y x if id == 4,    lpattern(dash)) ///
       (lfit ym xm, lpattern(solid)) ///
      , legend(order(2 "A" 3 "B" 4 "C" 5 "D" 1 "Means") size(small) title(Region:, size(small)) row(1)) ///
        xtitle("Component 1") ytitle("Component 2") ///
        title("Case B:" "Regional orthogonality," "individual correlation") ///
        xlabel(none) ylabel(none) name(caseb, replace)


// Case C
clear
set seed 1
set obs 3

generate  id = _n
generate u_i = 300 if id == 1
replace u_i = 0 if id == 2
replace u_i = -300 if id == 3

expand 150

bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70)                // Generate level-1 error term

generate y = x + u_i + e_ij
drop if y < 0 | x < 0

bys id: egen ym = mean(y)
bys id: egen xm = mean(x)

twoway (scatter y  x  if id == 1, msymbol(oh)) ///
       (scatter y  x  if id == 2, msymbol(dh)) ///
       (scatter y  x  if id == 3, msymbol(sh)) ///
       (scatter ym xm, msymbol(S)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
       (lfit ym xm, lpattern(solid)) ///
       , /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/ legend(off) ///
       xtitle("Component 1") ytitle("Component 2") ///
       title("Case C:" "Negative regional correlation," "positive individual correlation") ///
       xlabel(none) ylabel(none) ///
       name(casec, replace)

grc1leg casea caseb casec, legendfrom(caseb) span pos(6) row(1) name(combined, replace)