Reference Mummolo, Jonathan, and Erik Peterson. 2018. "Improving the Interpretation of Fixed Effects Regression Results." Political Science Research and Methods 6(4):829-835. doi: 10.1017/psrm.2017.44
set scheme lean1 webuse nlswork, clear xtset idcode year // Model of interest: xtreg ln_w ttl_exp age c.age#c.age tenure c.tenure#c.tenure not_smsa south i.year, fe local b : display %4.2f _b[ttl_exp] // Isolate relevant variation in the treatment: reghdfe ttl_exp, absorb(idcode year) residuals(relevantvariation) // Identify a plausible counterfactual shift in X given the data: center ttl_exp sum relevantvariation local sd1 : display %4.2f r(sd) sum ttl_exp local sd2 : display %4.2f r(sd) twoway (kdensity c_ttl_exp) /// (kdensity relevantvariation), /// ytitle("Density") xtitle(Treatment (centered around 0)) /// title("{bf:A} Treatment distribution", span bexpand justification(left)) /// legend(order(1 "Original distribution, {it:SD} = `sd2'" /// 2 "Distribution after FE, {it:SD} = `sd1'" ) /// pos(2) ring(0)) /// name(figure1A, replace) // Calculate range within individuals bysort idcode: egen max = max(ttl_exp) bysort idcode: egen min = min(ttl_exp) generate range = max - min drop max min qui sum range, detail local mean : display %4.2f r(mean) local p95 : display %4.2f r(p95) twoway (histogram range, freq) /// (scatteri 2000 `mean' (3) "Mean = `mean'", msymbol(i)) /// (scatteri 2000 `p95' (3) "95th percentile = `p95'" , msymbol(i)), /// title("{bf:B} Within-individual range of treatment", span bexpand justification(left)) /// ylabel(, format(%7.0gc)) ytitle(Frequency) /// xline(`mean') xline(`p95') /// xtitle(Within-individual range) legend(off) /// name(figure1B, replace) graph combine figure1A figure1B, col(2) // Calculate % of cases where treatment does not vary egen pickone = tag(idcode) fre range if pickone // In 12% of cases the treatment does not vary // Multiply the estimated coefficient of interest by the revised standard deviation display "Estimated coefficient of interest = `b'" display "Within-standard deviation = `sd1'" display "Product = " _skip(1) `sd1' * `b'
Showing posts with label twoway scatteri. Show all posts
Showing posts with label twoway scatteri. Show all posts
May 8, 2024
Mummolo & Peterson (2018): Fixed-effects regression
Jul 24, 2018
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)
Jul 22, 2018
Random graphs (134): Chi-squared function
twoway (function y = (x^(1 -1)*exp(-1*x/2))/((2^(1.0))*exp(lngamma(1.0))), range(0 20)) ///
(function y = (x^(1.5-1)*exp(-1*x/2))/((2^(1.5))*exp(lngamma(1.5))), range(0 20)) ///
(function y = (x^(2.5-1)*exp(-1*x/2))/((2^(2.5))*exp(lngamma(2.5))), range(0 20)) ///
(function y = (x^(5-1)*exp(-1*x/2))/((2^(5))*exp(lngamma(5))), range(0 20)) ///
(scatteri .5 0.1 "{it:df} = 2", msymbol(i) mlabpos(3)) ///
(scatteri .24 1.7 "{it:df} = 3", msymbol(i) mlabpos(3)) ///
(scatteri .15 5.0 "{it:df} = 5", msymbol(i) mlabpos(3)) ///
(scatteri .10 10.0 "{it:df} = 10", msymbol(i) mlabpos(3)) ///
, xtitle(Chi-squared) ytitle("") ylabel(, format(%6.1f)) legend(off) name(figu, replace) //ysize(6) xsize(6)
// Also search chidemo
Jul 12, 2018
Random graphs (133): Line plot
// Get divorce rate
sdmxuse data ESTAT, dataset(demo_ndivind) dimensions(A.DIVMARPCT.EU28) clear
drop freq
rename value divorce
tempname divorce
save `divorce', replace
// Get out-of-wedlock births
sdmxuse data ESTAT, dataset(demo_find) dimensions(A.NMARPCT.EFTA) clear
rename value emb
// Merge both data sets
merge 1:1 time using `divorce', nogenerate
destring time, replace
// Standardize
generate emb100 = (emb/5.6) * 100 // EMB in 1970
generate divorce100 = (divorce/11.3) * 100 // Divorce in 1970
// Plot
drop if time < 1970
twoway (line emb100 time) ///
(line divorce100 time) ///
(scatteri 410 2011 "Increase in divorce rate in Europe", msymbol(i) mlabpos(9)) ///
(scatteri 679 2016 "Increase in extramarital birth rate in Europe", msymbol(i) mlabpos(9)) ///
, ytitle("1970 = 100") xlabel(1970 (10) 2010) ylabel(100 (100) 700) xtitle("") legend(off) xsize(7)
Jan 7, 2018
Random graphs (124): Line plots
clear
input str100 grouping t1 t2
"Department managers" .22 .399
"Less than year" .35 .21
"Remaining non-managers" .31 .36
end
reshape long t, i(grouping) j(treatment)
twoway (connect t treatment if grouping == "Department managers", msymbol(o)) ///
(connect t treatment if grouping == "Less than year", msymbol(o)) ///
(connect t treatment if grouping == "Remaining non-managers", msymbol(o)) ///
(scatteri .399 2 "Department managers" ///
.21 2 "Less than year" ///
.36 2 "Remaining non-managers", msymbol(none)) ///
, legend(off) xtitle("Minimum wage raise") xlabel(0 " " 1 "Before" 2 "After" 3 " ", notick) ///
ytitle("Job satisfaction") ylabel(, format(%6.2f))
Jun 21, 2017
Chapter 4 of Singer and Willett's (2003) book on longitudinal data analysis
There is another take on the chapter here, but I like mine better.
use alcohol1_pp.dta, clear // Figure 4.1twoway (scatter alcuse age) (lfit alcuse age) /// if inlist(id, 4, 14, 23, 32, 41, 56, 65, 82), /// by(id, row(2) legend(off) note("")) /// xlabel(13 (1) 17) ytitle({it:ALCUSE}) xtitle({it:AGE}) /// name(figure41, replace) // Figure 4.2preserve bysort id: sample 32, count qui sum peer generate hipeer = (peer >= r(mean)) qui regress alcuse i.id##c.age predict predicted_alcuse label var predicted_alcuse "Predicted {it:ALCUSE}" xtline predicted_alcuse if coa == 0, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("{it:COA} = 0") ylabel(-1 (1) 4) /// name(coa0, replace) nodraw xtline predicted_alcuse if coa == 1, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("{it:COA} = 1") ylabel(-1 (1) 4) /// name(coa1, replace) nodraw xtline predicted_alcuse if hipeer == 0, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("Low {it:PEER}") ylabel(-1 (1) 4) /// name(hipeer0, replace) nodraw xtline predicted_alcuse if hipeer == 1, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("High {it:PEER}") ylabel(-1 (1) 4) /// name(hipeer1, replace) nodraw graph combine coa0 coa1 hipeer0 hipeer1, col(2) ysize(8) name(figure42, replace) restore // Table 4.1generate time = age - 14 // Create program that gets all the random slope parameters for esttab capture program drop randomslopetable program define randomslopetable eststo `1' qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance exp([lns1_1_1]_b[_cons]) * /// exp([lns1_1_2]_b[_cons]) qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance capture drop p1 // For R-squared based on observed and predicted outcomes qui predict p1 // For R-squared based on observed and predicted outcomes qui regress alcuse p1 // For R-squared based on observed and predicted outcomes local orsq = e(r2) // For R-squared based on observed and predicted outcomes estimates drop . // For R-squared based on observed and predicted outcomes estimates restore `1' // For R-squared based on observed and predicted outcomes qui estadd scalar orsq = `orsq' // For R-squared based on observed and predicted outcomes scalar v_e = exp(2*[lnsig_e]_b[_cons]) // For R-squared based on residual variance reduction scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar rsqe = (v_eumm - v_e) / v_eumm // For R-squared based on residual variance reduction of unconditional means model if !missing(`v2_ugm') { qui estadd scalar rsqv2 = (v2_ugm - v2) / v2_ugm // For R-squared based on residual variance reduction of unconditional means model qui estadd scalar rsqv1 = (v1_ugm - v1) / v1_ugm // For R-squared based on residual variance reduction of unconditional means model } end // This was based on: https://www.statalist.org/forums/forum/general-stata-discussion/general/1309801-saving-commands-as-macro-variables-local-and-global-to-make-do-files-shorter eststo clear mixed alcuse || id: , variance mle // Model A, unconditional means model eststo modela qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance scalar v_eumm = exp(2*[lnsig_e]_b[_cons]) // Residual variance for R-squared mixed alcuse time || id: time, cov(un) variance mle // Model B, unconditional growth model randomslopetable modelb scalar v2_ugm = exp(2*[lns1_1_1]_b[_cons]) // Slope variance scalar v1_ugm = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance mixed alcuse i.coa##c.time || id: time, cov(un) variance mle // Model C randomslopetable modelc mixed alcuse i.coa##c.time c.peer##c.time || id: time, cov(un) variance mle // Model D randomslopetable modeld mixed alcuse i.coa c.peer##c.time || id: time, cov(un) variance mle // Model E randomslopetable modele mixed alcuse i.coa c.cpeer##c.time || id: time, cov(un) variance mle // Model F randomslopetable modelf mixed alcuse ccoa c.cpeer##c.time || id: time, cov(un) variance mle // Model G randomslopetable modelg esttab /// , b(3) se(3) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define starts rename(ccoa 1.coa cpeer peer c.cpeer#c.time c.peer#c.time) /// // Align coefficients coeflabels(1.coa "COA" 1.coa#c.time "COA x time" peer "PEER" /// c.peer#c.time "PEER x time" _cons "Intercept" time "Time") /// // Label coefficients order(_cons 1.coa peer ) /// // Order coefficients stats(v_e v1 v2 cov orsq rsqe rsqv1 rsqv2 dev aic bic nc N, /// Add variance components to table fmt(3 3 3 3 3 3 3 3 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Cov(Init., Change)" /// "R-squared obs./pred." /// "R-squared var(Residual)" /// "R-squared var(Initial)" /// "R-squared var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C" "Model D" "Model E" "Model F" "Model G") /// keep(alcuse:) // Drop variance components in weird shapes // ICC of p. 96: estimates restore modela estat icc // Figure 4.3estimates restore modelb qui margins, at(time = (0 1 2)) marginsplot, xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") ylabel(0 .5 1 1.5 2) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Unconditional growth model", span) name(modelb, replace) nodraw estimates restore modelc qui margins, at(time = (0 1 2) coa = (0 1)) marginsplot, xlabel(-1 "13" 0 "14" .5 " " 1 "15" 1.5 " " 2 "16" 3 "17", ) ylabel(0 (1) 2, ) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Uncontrolled effects of {it:COA}", span) legend(off) /// addplot(scatteri .9 2 "{it:COA} = 0" 1.5 2 "{it:COA} = 1", msymbol(none)) /// name(modelc, replace) nodraw estimates restore modele margins, at(time = (0 1 2) coa = (0 1) peer = (.655 1.381)) marginsplot, xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17", format(%6.0f)) ylabel(0 1 2, ) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Controlled effects of {it:COA}", span bexpand) legend(off) /// addplot(scatteri .95 2 "{it:COA} = 0" /// 1.5 2 "{it:COA} = 1" /// .14 -0.85 "Low {it:PEER}" /// .64 -0.85 "High {it:PEER}" /// .72 -0.85 "Low {it:PEER}" /// 1.21 -0.85 "High {it:PEER}", msymbol(none)) /// name(modele, replace) nodraw graph combine modelb modelc modele, ycommon xcommon col(3) xsize(12) ysize(6) altshrink name(figure43, replace) // Figure 4.4preserve statsby, by(id) saving(temp, replace): regress alcuse time merge m:1 id using temp qui cor _b_cons coa local r = round(r(rho), .01) twoway (scatter _b_cons coa) /// (scatteri 4 .5 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(-.5 " " 0 1 1.5 " ") ytitle(Intercept) /// xtitle("{it:COA}") legend(off) name(g1, replace) nodraw qui cor _b_cons peer local r = round(r(rho), .01) twoway (scatter _b_cons peer) /// (scatteri 4 1 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(0 1 2 3) ytitle(Intercept) /// xtitle("{it:PEER}") legend(off) name(g2, replace) nodraw qui cor _b_time coa local r = round(r(rho), .01) twoway (scatter _b_time coa) /// (scatteri 4 .5 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(-.5 " " 0 1 1.5 " ") ytitle(Change) /// xtitle("{it:COA}") legend(off) name(g3, replace) nodraw qui cor _b_time peer local r = round(r(rho), .01) twoway (scatter _b_time peer) /// (scatteri 4 1 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(0 1 2 3) ytitle(Change) /// xtitle("{it:PEER}") legend(off) name(g4, replace) nodraw graph combine g1 g2 g3 g4, name(figure44, replace) restore // Figure 4.5estimates restore modelf predict e, resid predict e_time e_cons, reffects relevel(id) qnorm e, yline(0) ytitle(Var(residual)) name(g1, replace) nodraw qnorm e_cons, yline(0) ytitle(Var(intercept)) name(g2, replace) nodraw qnorm e_time, yline(0) ytitle(Var(change)) name(g3, replace) nodraw egen ze = std(e) egen ze_time = std(e_time) egen ze_cons = std(e_cons) scatter ze id, yline(0) ytitle("Standardized" "var(residual)") name(g4, replace) nodraw scatter ze_cons id, yline(0) ytitle("Standardized" "var(intercept)") name(g5, replace) nodraw scatter ze_time id, yline(0) ytitle("Standardized" "var(change)") name(g6, replace) nodraw graph combine g1 g2 g3 g4 g5 g6, cols(2) colfirst name(figure45, replace) ysize(8) // Figure 4.6twoway (scatter e time), ylabel(-2(1)2) ytitle("Var(residual)") /// xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") /// xtitle("{it:TIME}") /// yline(0) name(g1, replace) nodraw twoway (scatter e_cons coa), ylabel(-2(1)2) ytitle("Var(intercept)") /// xlabel(-.5 " " 0 1 1.5 " ") xtitle("{it:COA}") /// yline(0) name(g2, replace) nodraw twoway (scatter e_cons peer), ylabel(-2(1)2) ytitle("Var(intercept)") /// xlabel(0 1 2 3) xtitle("{it:PEER}") /// yline(0) name(g3, replace) nodraw twoway (scatter e_time coa), ylabel(-2(1)2) ytitle("Var(change)") /// xlabel(-.5 " " 0 1 1.5 " ") xtitle("{it:COA}") /// yline(0) name(g4, replace) nodraw twoway (scatter e_time peer), ylabel(-2(1)2) ytitle("Var(change)") /// xlabel(0 1 2 3) xtitle("{it:PEER}") /// yline(0) name(g5, replace) nodraw graph combine g1 g2 g3 g4 g5, cols(2) hole(2) name(figure46, replace) ysize(8) // Figure 4.7qui reg alcuse time coa cpeer c.cpeer#c.time predict pa estimates restore modelf gen bayes0 = _b[_cons] + _b[1.coa]*coa + _b[cpeer]*cpeer + e_cons gen bayes1 = _b[time] + _b[c.cpeer#c.time]*cpeer + e_time gen bayes = bayes0 + bayes1 * time twoway (scatter alcuse time) /// (lfit alcuse time) /// (line pa time, sort) /// (line bayes time, sort) /// if inlist(id, 4, 14, 23, 32, 41, 56, 65, 82), by(id, col(4) note("")) /// ytitle({it:ALCUSE}) xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") /// legend(order(1 "Data points" 2 "OLS" /// 3 "Population average" 4 "Bayes estimate") col(4)) /// xtitle("{it:TIME}") name(figure47, replace)
Reference
Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis. Modeling Change and Event Occurrence. Oxford University Press. doi: 10.1093/acprof:oso/9780195152968.001.0001
Jun 20, 2017
Random graphs (100): Labeling lines directly
import delimited une_rt_a_1_Data.csv, clear
/*
DATASET: Unemployment by sex and age - annual average [une_rt_a]
LAST UPDATE: 14.06.17 13:10:19
EXTRACTION DATE: 18.06.17 11:35:33
SOURCE OF DATA: Eurostat
*/
drop sex age unit flagandfootnotes
// Generate country variable
kountryadd "Germany (until 1990 former territory of the FRG)" to "Germany" add
kountry geo, from(other) stuck marker
ren _ISO3N_ country
kountry country, from(iso3n) to(iso2c)
ren _ISO2C_ country_str
replace country_str = "UK" if country_str == "GB"
list country_str geo
// Fix unemployment rate
replace value = "" if value == ":" // Fix missing data indicator
destring value, replace // Convert to numeric
twoway (line value time if country_str == "UK") ///
(line value time if country_str == "US") ///
(scatteri 4 2016 "United States" 3.5 2016 "United Kingdom", msymbol(none)) ///
, legend(off) ///
xtitle("") ytitle("Male unemployment rate, 25-74 y.") ///
xlabel(1985(5)2015) ///
xscale(range(1983 2023)) ///
note(" " "{it: Source:} Eurostat, une_rt_a, date of extraction: 2017-06-18", span)
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)
Jun 18, 2017
Random graphs (98): Shaded areas
import delimited une_rt_a_1_Data.csv, clear
/*
DATASET: Unemployment by sex and age - annual average [une_rt_a]
LAST UPDATE: 14.06.17 13:10:19
EXTRACTION DATE: 18.06.17 11:35:33
SOURCE OF DATA: Eurostat
*/
drop sex age unit flagandfootnotes
// Generate country variable
kountryadd "Germany (until 1990 former territory of the FRG)" to "Germany" add
kountry geo, from(other) stuck marker
ren _ISO3N_ country
kountry country, from(iso3n) to(iso2c)
ren _ISO2C_ country_str
replace country_str = "UK" if country_str == "GB"
list country_str geo
// Fix unemployment rate
replace value = "" if value == ":" // Fix missing data indicator
destring value, replace // Convert to numeric
// Plot
twoway (scatteri 12 2008 12 2012, recast(area) bcolor(gs14)) ///
(line value time if country_str == "UK") ///
(line value time if country_str == "DE") ///
(line value time if country_str == "US") ///
, legend(order(2 "UK" 3 "Germany" 4 "US") ///
pos(7) ring(0)) ///
xtitle("") ytitle("Male unemployment rate, 25-74 y.") xtick(1983(1)2016) ///
xlabel(1985(5)2015) ///
note(" " "{it: Source:} Eurostat, une_rt_a, date of extraction: 2017-06-18", span)
Subscribe to:
Posts (Atom)