* HINT: * The program assumes that all files are stored in * C:\monographregression\computercode * change directory if code is located elsewhere * to store the figures a subfolder `figures' has to be * created clear cd C:\monographregression\computercode set scheme s1mono * PART I data preparation (effect coding, polynomials etc.) use "data\stata\beach.dta" * Drop observations where the pH-value is missing drop if ph==. * Create orthogonal polynomials of year and of age orthpoly year , generate(yearo yearo2 yearo3) degree(3) orthpoly age , generate(ageo ageo2 ageo3) degree(3) * Generate a categorical variable for the degree of defoliation generate defol3 = 1 replace defol3 = 2 if defol >= 12.5 replace defol3 = 3 if defol >= 50 * Generate effect-coded variables for the thickness of humus layer generate humus0 = 0 replace humus0 = 1 if humus==0 replace humus0 = -1 if humus==1 generate humus2 = 0 replace humus2 = 1 if humus==2 replace humus2 = -1 if humus==1 generate humus3 = 0 replace humus3 = 1 if humus==3 replace humus3 = -1 if humus==1 generate humus4 = 0 replace humus4 = 1 if humus>3 replace humus4 = -1 if humus==1 * Generate effect-coded variables for the level of soil moisture generate watermoisture1 = 0 replace watermoisture1 = 1 if watermoisture==1 replace watermoisture1 = -1 if watermoisture==2 generate watermoisture3 = 0 replace watermoisture3 = 1 if watermoisture==3 replace watermoisture3 = -1 if watermoisture==2 * Generate effect-coded variables for fraction of alkali-ions in soil generate alkali1 = 0 replace alkali1 = 1 if alkali==1 replace alkali1 = -1 if alkali==2 generate alkali3 = 0 replace alkali3 = 1 if alkali==3 replace alkali3 = -1 if alkali==2 generate alkali4 = 0 replace alkali4 = 1 if alkali==4 replace alkali4 = -1 if alkali==2 * Generate an effect-coded variable for the fertilization replace fert = -1 if fert==0 * Generate an effect-coded variable for the type of forest replace type = -1 if type==0 * PART II models * Cumulative logit model ologit defol3 ageo ageo2 ageo3 yearo yearo2 yearo3 /* */ gradient canopyd alt depth ph type fert humus0 humus2 humus3 humus4 /* */ watermoisture1 watermoisture3 alkali1 alkali3 alkali4 * Save the estimated coefficients matrix b = e(b) * Calculate the estimated effects of age and year generate fage = b[1,1]*ageo + b[1,2]*ageo2 + b[1,3]*ageo3 generate fyear = b[1,4]*yearo + b[1,5]*yearo2 + b[1,6]*yearo3 * Plot of the effects of age and year graph twoway line fage age, sort clpattern(solid) lwidth(medthick) xsize(4.5) ytitle("") xtitle("age", size(medium)) /* */ ylab(-4(1)2, labsize(medium) nogrid) xlab(5(25)230, labsize(medium)) title("") legend(off) subtitle("Effect: age") graph export figures\fage_param.eps, replace graph twoway line fyear year, sort clpattern(solid) lwidth(medthick) xsize(4.5) ytitle("") xtitle("year", size(medium)) /* */ ylab(-0.75(0.25)0.5 , labsize(medium) nogrid) xlab(1983(7)2004, labsize(medium)) title("") legend(off) subtitle("Effect: calendar time") graph export figures\fyear_param.eps, replace * Duplicate the data sort age expand 2 * Set orthogonal polynomials of year to zero replace yearo = 0 if _n >= 1794 replace yearo2 = 0 if _n >= 1794 replace yearo3 = 0 if _n >= 1794 * Set the continuous covariates to their mean values sum gradient if _n < 1794 replace gradient = _result(3) if _n >= 1794 sum canopyd if _n < 1794 replace canopyd = _result(3) if _n >= 1794 sum alt if _n < 1794 replace alt = _result(3) if _n >= 1794 sum depth if _n < 1794 replace depth = _result(3) if _n >= 1794 sum ph if _n < 1794 replace ph = _result(3) if _n >= 1794 * Set the categorical covariates to their mode category replace type = 1 if _n >= 1794 replace fert = -1 if _n >= 1794 replace humus0 = -1 if _n >= 1794 replace humus2 = -1 if _n >= 1794 replace humus3 = -1 if _n >= 1794 replace humus4 = -1 if _n >= 1794 replace watermoisture1 = -1 if _n >= 1794 replace watermoisture3 = -1 if _n >= 1794 replace alkali1 = -1 if _n >= 1794 replace alkali3 = -1 if _n >= 1794 replace alkali4 = -1 if _n >= 1794 * Cumulative logit model ologit defol3 ageo ageo2 ageo3 yearo yearo2 yearo3 /* */ gradient canopyd alt depth ph type fert humus0 humus2 humus3 humus4 /* */ watermoisture1 watermoisture3 alkali1 alkali3 alkali4 if _n < 1794 * Save the estimated probabilities predict mu1 mu2 mu3 * Plot the estimated probabilities scatter mu1 age if _n>=1794 , sort clpattern(l) m(i) c(l) lwidth(medthick) xsize(4.5) ytitle("") xtitle("age", size(medium)) /* */ ylab(0.4(0.1)1, labsize(medium) nogrid) xlab(5(25)230, labsize(medium)) title("") legend(off) subtitle("Estimated probability for no damage") graph export figures\fage_mu1.eps, replace scatter mu2 age if _n>=1794 , sort clpattern(l) m(i) c(l) lwidth(medthick) xsize(4.5) ytitle("") xtitle("age", size(medium)) /* */ ylab(0(0.1)0.6, labsize(medium) nogrid) xlab(5(25)230, labsize(medium)) title("") legend(off) subtitle("Estimated probability for weak damage") graph export figures\fage_mu2.eps, replace scatter mu3 age if _n>=1794 , sort clpattern(l) m(i) c(l) lwidth(medthick) xsize(4.5) ytitle("") xtitle("age", size(medium)) /* */ ylab(, labsize(medium) nogrid) xlab(5(25)230, labsize(medium)) title("") legend(off) subtitle("Estimated probability for severe damage") graph export figures\fage_mu3.eps, replace