Showing posts with label twoway lfit. Show all posts
Showing posts with label twoway lfit. Show all posts

Jul 24, 2018

Random graphs (137): Logistic regression

clear

// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly
qui sum questions, detail
generate pass = (questions >= r(p75)) // Passing the exam

// 1) Histogram of outcome variable
twoway (histogram pass, discrete percent), ///
        xlabel(0 "[0] Failed" 1 "[1] Passed") xtitle("") ///
        ytitle("Percent of students") xsize(4) ysize(4) name(figure5, replace)

// 2) Scatterplot
twoway (scatter pass hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure6, replace)    

// 3) Scatterplot with regression line
regress pass hours
local intercept = round(_b[_cons], .01)
local x = round(_b[hours], .01)      
       
twoway (scatter pass hours) ///
       (lfit pass hours, lpattern(solid) range(1 5)), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(.8 2 "y = `intercept' + `x' x + e", size(large)) ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure7, replace)

// 4) Logit curve
logit pass hours
predict yhat   

twoway (scatter pass hours) ///
       (line yhat hours, lpattern(solid) sort), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure8, replace)

graph combine figure5 figure6 figure7 figure8, ///
              col(2) xsize(8) ysize(8) altshrink name(figures58, replace)

Random graphs (135): Scatterplot with OLS regression

clear

// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly

// 1) Basic scatterplot
twoway (scatter questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Number of questions answered correctly") ///
        name(figure1, replace) ylabel(0 (5) 15) xsize(4) ysize(4)

// 2) Scatterplot with regression line
twoway (scatter questions hours) ///
       (lfit questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 (5) 15) legend(off) name(figure2, replace) xsize(4) ysize(4)

// 3) Scatterplot with regression line and equation
regress questions hours
local intercept = round(_b[_cons], .1)
local x = round(_b[hours], .1)

twoway (scatter questions hours) ///
       (lfit questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(14 2 "y = `intercept' + `x' + e", size(large)) ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 (5) 15) legend(off) name(figure3, replace) xsize(4) ysize(4)

// 4) Scatterplot with regression line and equation and labels for components
twoway (scatter questions hours) ///
       (lfit questions hours, lpattern(solid) range(0 5)) ///
       (function y = 1.6 + 2.5, range(1 2)) ///
       (function y = 2, range(4.1 6.6) horizontal lpattern(solid) lcolor(red)) ///
       (scatteri 1.6 0 (3) "Intercept", msymbol(o) mlabcolor(red)) ///
       (scatteri 5 2 "Slope", msymbol(i) mlabcolor(red)), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(14 2 "y = `intercept' + `x' x + e", size(large)) ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 1.6 5 10 15) legend(off) name(figure4, replace) xsize(4) ysize(4)
    
graph combine figure1 figure2 figure3 figure4, col(2) xsize(8) ysize(8) altshrink name(figures14, replace)

Jun 19, 2017

Random graphs (99): Scatterplots

regress treatments index
local r2 = round(e(r2), .01) // Round coeff.
twoway (scatter treatments index) ///
       (lfit treatments index) ///
       (scatter treatments index if inlist(cntry, "IT", "AT", "UA", "BE", "DK", "RU"), msymbol(o)) ///
      , xtitle(ART comprehensive availability) ytitle("ART treatments" "per million women 15-44 y.") ///
        title("{bf:B}", justification(left) bexpand span) ///
        xlabel(0 (1) 9) ///
        text(1500 6.6 "Denmark") ///   
        text(1500 8.6 "Belgium") ///   
        text(470  0.4 "Italy") ///
        text(400  1.6 "Austria") ///
        text(120  7.5 "Ukraine") ///
        text(210  8.0  "Kazakhstan", placement(west)) ///
        text(150  8.0 "Russia", placement(east)) ///
     legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(7)) name(policy, replace)

Dec 17, 2014

Random graphs (39): Plotting on a log axis


Rather than looking at a variable in absolute terms, namely GDP per capita in the upper panel of the Figure, it is also sometimes helpful to look at it in terms of percentage increases, as in the lower panel. Equal distances on the x-axis refer to equal percentage increases in GDP per capita. Four equidistant points on the x-axis are labeled, each indicating a fourfold increase in GDP.

use wvs2005_v20090901a.dta, clear

// Country variable
// 1) Turn into string
decode v2, gen(ctry)

// 2) Abbreviate
kountry ctry, from(other) stuck marker
ren _ISO3N_ cntry
kountry cntry, from(iso3n) to(iso2c)
ren _ISO2C_ country

// Generate outcome: % in good health
generate goodhealth = 100 if v11 <= 2
replace  goodhealth =   0 if v11  > 2
replace  goodhealth =   . if v11 == .

// Collapse data set
collapse (mean) goodhealth [pw = v259], by(country)

// Generate year variable
gen year = 2006

// Preserve collapsed data set
preserve

// Get GDP from World Bank data base
wbopendata, language(en - English) country() topics() indicator(NY.GDP.PCAP.PP.CD) clear long

// Fix obtained data set
ren ny_gdp_pcap_pp_cd gdp
ren iso2code country
keep if year == 2006
keep gdp country
drop if country == ""

// Save obtained GDP data
tempfile gdp
save `gdp'

// Restore 
restore

// Merge GDP with collapsed data set
merge m:1 country using `gdp'
keep if _merge == 3
drop _merge

// Create labels with thousand separator
label define gdp 20000 "20,000" ///
                 40000 "40,000" ///
                 60000 "60,000"
label val gdp gdp

// Plot on unlogged axis
twoway (scatter goodhealth gdp, mlabel(country) mlabpos(0) msymbol(none)) ///
       (lfit goodhealth gdp) ///
      , legend(off) xtitle("GDP per capita, 2006, PPP in current international $") ///    
        ytitle("% in good health") ///
        xlabel(0(20000)60000, valuelabels) ///
        name(unlogged, replace)

// Generate logged variable
gen loggdp = log(gdp) * 1000  // Multiply by 1,000 because only integers can be labeled

// Generate numbers for labeling
//  Round them to three decimal digits, then multiply by 1,000 to get integers
local log1 = round(log(1000), .001) * 1000
local log2 = round(log(1000 * 4), .001) * 1000
local log3 = round(log(1000 * 4 * 4), .001) * 1000
local log4 = round(log(1000 * 4 * 4 * 4), .001) *1000

*di `log1' _skip(2) `log2' _skip(2) `log3' _skip(2) `log4' 

// Create labels for numbers to be labeled
label define loggdp `log1' "1,000" ///
                    `log2' "4,000" ///
                    `log3' "16,000" ///                    
                    `log4' "64,000"
label value loggdp loggdp

// Plot on log axis
twoway (scatter goodhealth loggdp, mlabel(country) mlabpos(0) msymbol(none)) ///
       (lfit goodhealth loggdp) ///
      , legend(off) xtitle("GDP per capita, 2006, PPP in current international $") ///   
        ytitle("% in good health") ///
        xlabel(`log1' `log2' `log3' `log4', valuelabels) ///
        name(logged, replace)

// Combine plots
graph combine unlogged logged, col(1) ysize(8)

Oct 24, 2014

Random graphs (34): Country scatterplots

use ART.dta, replace

// Fix missing values
recode avgemtrans (-1 -3 = .)
recode emhum      (-1 = .)

// Calculate r-squared
qui regress avgemtrans emhum
local r2 = round(e(r2), .01) // Round coeff.

// Create first panel
twoway (scatter avgemtrans emhum, mlab(countrycode) mlabpos(0) msymbol(i)) ///
       (lfit avgemtrans emhum) ///
    , xlabel(, format(%6.1f)) ylabel(, format(%6.1f)) ///
    ytitle("Avg. no. of fresh non-donor embryos" "per fresh embryo transfer cycle") ///
    xtitle("Avg. level disagreement that" "embryo is a human being") ///
       legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
    name(figure2a, replace)

// Fix missing values
recode afford clonehelpinf (-1 = .)

// Listwise deletion to control y-axis range
preserve
keep if !missing(afford, clonehelpinf)

// Calculate r-squared    
regress afford clonehelpinf
local r2 = round(e(r2), .01) // Round coeff.

// Create second panel
twoway (scatter afford clonehelpinf, mlab(countrycode) mlabpos(0) msymbol(i)) ///
       (lfit afford clonehelpinf) ///
    , xlabel(, format(%6.1f)) ylabel(0 (5) 25, format(%6.1f)) ///
    ytitle("Affordability (net cost of a fresh ART cycle" "as % of annual disposable income)") ///
    xtitle("Avg. level disagreement with" "cloning to help infertile couples") ///
       legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
    name(figure2b, replace)
restore

// Create final Figure
graph combine figure2a figure2b, row(1)

Jul 1, 2014

Random graphs (23): Between and within regression

clear

input j i xij x_j yij y_j
// Data from Snijders and Bosker (1999, Table 3.2):
1 1 1 2 5 6
1 2 3 2 7 6
2 1 2 3 4 5
2 2 4 3 6 5 
3 1 3 4 3 4
3 2 5 4 5 4
4 1 4 5 2 3
4 2 6 5 4 3
5 1 5 6 1 2
5 2 7 6 3 2
end


// Plot as Figure 3.4 of Snijders ans Bosker (1999):
twoway (scatter yij xij) ///
       (lfit yij xij) ///
       (lfit y_j x_j) /// 
       (lfit yij xij if j == 1, lpattern(dash_dot)) ///
       (lfit yij xij if j == 2, lpattern(dash_dot)) ///
       (lfit yij xij if j == 3, lpattern(dash_dot)) ///
       (lfit yij xij if j == 4, lpattern(dash_dot)) ///
       (lfit yij xij if j == 5, lpattern(dash_dot)) ///
     , legend(label(2 "Total regression") ///
       label(3 "Between regression") ///
       label(4 "Within regression") ///
       order(2 3 4) ///
       pos(1) ring(0)) ///
       xtitle(X) ytitle(Y) ///
       title("Within, between, and total relations") /// 
       xlabel(none) ylabel(none) ///
       name(wbt_legend, replace)
    

// Lines directly labeled
twoway (scatter yij xij) ///
       (lfit yij xij) ///
       (lfit y_j x_j) ///
       (lfit yij xij if j == 1, lpattern(dash_dot)) ///
       (lfit yij xij if j == 2, lpattern(dash_dot)) ///
       (lfit yij xij if j == 3, lpattern(dash_dot)) ///
       (lfit yij xij if j == 4, lpattern(dash_dot)) ///
       (lfit yij xij if j == 5, lpattern(dash_dot)) ///
     , legend(off) ///
       text(1.8 6.5 "Between") ///
       text(3.0 7.5 "Total") ///
       text(7.0 3.5 "Within") ///
       text(5.0 5.5 "Within") ///
       xtitle(X) ytitle(Y) ///
       title("Within, between, and total relations") ///
       xlabel(none) ylabel(none) ///
       name(wbt_nolegend, replace)

graph combine wbt_legend wbt_nolegend, col(1) ysize(8) xcommon ycommon
// ysize(4) and xsize(5.5) are defaults

Reference

 Snijders, Tom, and Roel Boskers. 1999. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling. Sage.

Jul 9, 2013

Morgan and Winship's (2007) example for bias due to conditioning on a collider in Stata

Morgan and Winship (2007: p. 66) illustrate Pearl's (2009) concern about conditioning on a collider variable using a simple example.

The general problem of conditioning on a collider is as follows. Consider three variables A, B, and C, with both A and B being causes of C: A → C ← B. (Formally, any variable C that has two arrows pointing to it along a given path is a collider.) Unlike a confounder (an uncontrolled common cause of A and B), a collider does not induce a zero-order correlation between A and B. However, when handled inappropriately, a collider can induce a conditional correlation between A and B.

Morgan and Winship's example shows just that: A college admits applicants based on their SAT scores and ratings of their motivation based on an interview. Those in the top 15 per cent of the sum of SAT and motivation ratings are being admitted. SAT scores and motivation ratings are largely uncorrelated.

// Generate and label two variables
drawnorm sat motivation, ///
         n(250) ///
         means(.007, -.053) ///
         sds(1.01, 1.02) ///
         corr(1, .035, 1) cstorage(lower) ///
         clear seed(1)
label var sat        "SAT"
label var motivation "Motivation"

// Only the 15 per cent at the top are admitted
gen admission_sc = sat + motivation
_pctile admission_sc, percentiles(85)

gen admission = (admission_sc > r(r1))
label var admission "Admission status"
label define admission 1 "Admitted applicant" ///
                       0 "Rejected applicant"
label val admission admission
drop admission_sc

// Plot as in Morgan and Winship, p. 67:
twoway (scatter motivation sat if admission == 1) ///
       (scatter motivation sat if admission == 0) ///
       , ///
       legend(label(1 "Admitted applicants") ///
              label(2 "Rejected applicants") ///
              pos(5) ring(0)) ///
       ylabel(none) xlabel(none) ///
       name(collider1, replace)
// Enhanced plot with fitted lines and correlations
quietly cor motivation sat
local r_overall = round(r(rho), .01)
quietly cor motivation sat if admission == 1
local r_admitted = round(r(rho), .01)
quietly cor motivation sat if admission == 0
local r_rejected = round(r(rho), .01)
    
twoway (scatter motivation sat if admission == 1) ///
       (scatter motivation sat if admission == 0) ///
       (lfit motivation sat) ///
       (lfit motivation sat if admission == 1) ///
       (lfit motivation sat if admission == 0) ///
       , ///
       legend(label(1 "Admitted applicants") ///
              label(2 "Rejected applicants") ///
              label(3 "Overall fit, {it:r} = `r_overall'") ///
              label(4 "Fit for admitted, {it:r} = `r_admitted'") ///
              label(5 "Fit for rejected, {it:r} = `r_rejected'") ///
              pos(5) ring(0)) ///
       ylabel(none) xlabel(none) ytitle("Motivation")

What the Figures show is that the very small correlation between motivation and SAT score for the overall group turns out to be much larger when conditioning for admission status.

This also shows in an OLS regression:

regress motivation sat, beta
estimates store m1
regress motivation sat admission, beta
estimates store m2

estimates table m1 m2, b(%7.2f) se(%7.2f) stats(N) label

----------------------------------------------
                Variable |   m1        m2     
-------------------------+--------------------
                     SAT |    0.10     -0.21  
                         |    0.06      0.06  
        Admission status |              1.60  
                         |              0.17  
                Constant |   -0.12     -0.36  
                         |    0.06      0.06  
-------------------------+--------------------
                       N |     250       250  
----------------------------------------------
                                  legend: b/se

Cole et al. (2010) present additional illustrations for this problem.

References


Cole, Stephen R., Robert W. Platt, Enrique F. Schisterman, Haitao Chu, Daniel Westreich, David Richardson, and Charles Poole. 2010. "Illustrating Bias Due to Conditioning on a Collider." International Journal of Epidemiology 39(2):417-420. doi: 10.1093/ije/dyp334

Morgan, Stephen L., and Christopher Winship. 2007. Counterfactuals and Causal Inference. Methods and Principles for Social Research. Cambridge University Press.

Pearl, Judea. 2009. Causality. Models, Reasoning, and Inference, 2nd ed. Cambridge University Press.

Oct 5, 2012

Random graphs (4): Plotting intercept and slope variation


use V3 V51 ISCO88 using "C:\Users\User\work\data sets\issp - work orientations\work orientations iii (2005)\ZA4350_F1.dta", clear

// Generate variables
gen jobsatisf = 7 - V51
iskoisei isei, isko(ISCO88)
drop V51 ISCO88
preserve

// Calculating intercept and slope variation using OLS statsby inter = _b[_cons] /// slope = _b[isei] /// , by(V3) /// saving(ols, replace): /// regress jobsatisf isei merge m:1 V3 using ols drop _merge // Visualizing intercept and slope variation gen yhat_ols = inter + slope*isei separate jobsatisf, by(V3) separate yhat_ols, by(V3) twoway (line yhat_ols1-yhat_ols43 isei, sort(V3 isei)) /// (lfit jobsatisf isei, clwidth(vvthick) clcolor(black)) /// , legend(off) ytitle("Job satisfaction") xtitle("ISEI") /// xlabel(16 25 50 75 90) ylabel(,format(%6.1f)) /// caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) /// name(one, replace) restore
// Calculating intercept and slope variation using a random effects model center isei, inplace mixed jobsatisf isei || V3: isei, var predict u1 u0, reffects // Visualizing intercept and slope variation gen predRandomSlope = (_b[_cons] + u0) + ((_b[isei] + u1) * isei) twoway (line predRandomSlope isei, connect(ascending) sort(V3 isei)), /// ytitle("Job satisfaction") /// xtitle("ISEI (centered)") /// xlabel(-25 0 25 50) /// ylabel(,format(%6.1f)) /// caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) /// name(two, replace)