/* Paper: Effect of tobacco control policies on the Swedish Smoking Quitline using intervention time series analysis. BMJ Open. Authors: Xingwu Zhou, Alessio Crippa, Anna-Karin Danielsson, Rosaria Galanti, and Nicola Orsini Affiliation: Department of Public Health Sciences, Karolinska Institutet Goal: This do-file can be useful to reproduce the results Author: Orsini N Date: October 17, 2019 */ * Open the dataset preserve use http://www.stats4life.se/data/quitline, clear tsset time, monthly * Evaluation of the effect of the EU Directive on May 2016 (Figure 1) keep if inrange(time, 625, 691) * Generate transformations of time to be used in the Poisson model * Splines of degree 0 and 1 with a knot at the intervention time of 676 (May 2016) gen s0676 = (time>676) gen s1676 = (time>676)*(time-676) * Fourier transformations * degrees variable for time divided by the number of time points in a year (i.e. 12 for months) gen degrees=(time/12)*360 fourier degrees, n(3) * Poisson model considering seasonality, change in level at the intervention, and change in linear trends glm call cos_* sin_* time s0676 s1676, fam(poisson) lnoffset(pop) link(log) scale(x2) * Shift at the intervention lincom s0676, eform * Linear trend before intervention (per year) lincom time*12, eform * Linear trend after intervention (per year) lincom (time+s1676)*12, eform * Graph the predicted rates over a fine grid of time values set obs 500 range timef 625 691 capture drop cos_* sin_* capture drop degreesf gen degreesf=(timef/12)*360 fourier degreesf, n(3) capture drop s0676 s1676 gen s0676 = (timef>676) gen s1676 = (timef>676)*(timef-676) predictnl fitc3 = _b[_cons] + _b[cos_1]*cos_1+ _b[cos_2]*cos_2 + _b[cos_3]*cos_3 /// + _b[sin_1]*sin_1+ _b[sin_2]*sin_2 + _b[sin_3]*sin_3 + _b[time]*timef replace fitc3 = exp(fitc3)*10^5 predictnl fitp3 = _b[_cons] + _b[cos_1]*cos_1+ _b[cos_2]*cos_2 + _b[cos_3]*cos_3 /// + _b[sin_1]*sin_1+ _b[sin_2]*sin_2 + _b[sin_3]*sin_3 + _b[time]*timef + _b[s0676]*s0676+ _b[s1676]*s1676 replace fitp3 = exp(fitp3)*10^5 * Fix transformations at April -.0608104 -.9926042 .1815316 .9981493 -.1213956 -.9833851 predictnl fitc31 = _b[_cons] + _b[cos_1]* -.0608104 + _b[cos_2]* -.9926042 + _b[cos_3]* .1815316 /// + _b[sin_1]* .9981493 + _b[sin_2]*-.1213956 + _b[sin_3]* -.9833851 + _b[time]*timef predictnl fitc32 = _b[_cons] + _b[cos_1]* -.0608104 + _b[cos_2]* -.9926042 + _b[cos_3]* .1815316 /// + _b[sin_1]* .9981493 + _b[sin_2]*-.1213956 + _b[sin_3]* -.983385 + _b[time]*timef + _b[s0676]*s0676+ _b[s1676]*s1676 replace fitc31 = exp(fitc31)*10^5 replace fitc32 = exp(fitc32)*10^5 gen low = 40 gen high = 120 twoway /// (rarea low high timef if inrange(timef,676, .), color(gs14) ) /// (scatter rate time, ms(o) c(l) mc(black) msize(small) lp(dot) lw(vthin) lc(black)) /// (line fitp3 timef if timef >= 676.6, sort lc(blue) ) /// (line fitc3 timef if timef <= 676, sort lc(blue)) /// (line fitc31 timef if timef <= 676, sort lc(red)) /// (line fitc31 timef if timef >= 677, sort lp(dash) lc(red)) /// (line fitc32 timef if timef >= 676.8, sort lc(red)) /// , plotregion(style(none)) scheme(s1mono) /// ytitle("Calling rates per 100,000 smokers") /// xtitle("Calendar time") ylabel(#5, angle(horiz)) /// xlabel(625(3)691, angle(45) format(%tm) ) xmtick(625(1)691) legend(off) name(figure1, replace) restore