The empirical time trends shown in Figure 2.4 on page 8 indicated very dissimilar trends for the three categories. In contrast, the regression model (12.1) assumes a common time trend and, in fact, the estimates shown in Figure 12.1 seem to be dominated by the trend for slightly damaged trees. To relax the assumption of a common time trend, we replaced the predictor in (12.1) by
ηit(r)=θ(r)−h
f1(r)(t) +f2(ait) +fspat(si) +u0itγi
, (12.2)
where f1(r)(t) is a category-specific time trend. In contrast to (12.1), this model does not contain an interaction effect, since the category-specific time trend would also imply a category-specific interaction making model (12.2) hardly identifiable from the present data. We experimented with discrete interactions based on a categorized age effect but, due to the small number of cases in category ’3’, these models turned out to be either not
−2−1012
1983 1990 1997 2004
calendar time
time trend 1
−2−1012
1983 1990 1997 2004
calendar time
time trend 2
Figure 12.6: Forest health data: Posterior mode estimates of category-specific time trends when the spatial effect is modeled by a Markov random field.
identifiable or not interpretable. Therefore, we decided to exclude the interaction term and to estimate the main effects model (12.2) for illustration purposes. The spatial effect is again assumed to follow a Markov random field prior.
Within BayesX, category-specific effects are requested by adding the additional option catspecificto the respective effect. In our example, the specificationtime(psplinerw2) has to replaced by time(psplinerw2,catspecific). Note that parametric effects can also be specified to be category-specific, i. e. the specification of x(catspecific) corre-sponding to effects xγ(r) is also allowed.
Figure 12.6 visualizes the two estimated trend functions. The trend f1(1)(t), which influ-ences the first threshold, almost equals the overall trend obtained in model (12.1) except for the narrower credible intervals. In contrast, the second trend f1(2)(t) exhibits a dif-ferent shape. Starting from a negative value, it increases over most of the observation period before showing a slight decrease at the end of the nineties. This corresponds very well to the empirical trend of category ’3’ which shows a similar pattern (see Figure 2.4 on page 8).
13 Simulation studies for multicategorical responses
To investigate the performance of the presented mixed model based approach to categor-ical structured additive regression, we conducted several simulation studies based on a multinomial logit model and a cumulative probit model with three categories. In either model the predictors were defined to be the sum of a nonparametric effect and a spatial effect, see Figures 13.1 and 13.2 for detailed descriptions of the simulation design.
-1.7 0 1.4
Category 1
-1.7 0 1.4
Category 2
- Predictor:
ηi(r) =f1(r)(xi) +f2(r)(si) - Category 1:
f1(1)(x) = sin[π(2x−1)]
f2(1)(s) = −0.75|sx|(0.5 +sy) - Category 2:
f1(2)(x) = sin[2π(2x−1)]
f2(2)(s) = 0.5(sx+sy)
- xis chosen from an equidistant grid of 100 values between -1 and 1.
- (sx, sy) are the centroids of the 124 districts s of the two southern states of Germany (see Figures).
Figure 13.1: Simulation design for the multinomial logit model.
13.1 Comparison of different modeling approaches
The first aim of this simulation study was the comparison of different parameterizations of the spatial effect and of different approaches to the estimation of categorical STAR models. In total, 250 simulation runs with n= 500 observations were generated from the multinomial logit model described in Figure 13.1. We used cubic P-splines with second order random walk penalty and 20 knots to estimate effects of the continuous covariate.
The spatial effect was estimated either by a Markov random field, a (full) Gaussian ran-dom field or a two-dimensional P-spline (based on 10×10 inner knots). Concerning the competing fully Bayesian approach by Fahrmeir & Lang (2001b) and Brezger & Lang
-1.7 0 1.4
- Predictor:
ηi(r) =θ(r)−f1(xi)−f2(si) - Functions:
f1(x) = sin[π(2x−1)]
f2(s) = 0.5(sx+sy)
- xis chosen from an equidistant grid of 100 values between -1 and 1.
- (sx, sy) are the centroids of the 124 districts s of the two southern states of Germany (see Figure).
Figure 13.2: Simulation design for the cumulative probit model.
(2005), inverse Gamma priors IG(a, b) witha=b= 0.001 were assigned to the variances.
In the fully Bayesian framework, the GRF approach was computationally to demanding since it requires the inversion of a full precision matrix for the spatial effect in each it-eration. Therefore, we excluded the fully Bayesian GRF approach from the comparison.
As a further competitor, we utilized the R-implementation of the procedure polyclass de-scribed in Kooperberg et al. (1997), where nonparametric effects and interaction surfaces are modeled by linear splines and their tensor products. Smoothness of the estimated curves is not achieved by penalization but via stepwise inclusion and deletion of basis functions based on AIC.
The results of the simulation study can be summarized as follows:
• Generally, REML estimates have somewhat smaller median MSE than their fully Bayesian counterparts, with larger differences for spatial effects (see Figures 13.3a to 13.3d).
• Estimates for the effects of the continuous covariate are rather insensitive with respect to the model choice for the spatial effect (Figures 13.3a and 13.3b).
• Two-dimensional P-splines lead to the best fit for the spatial effect although data are provided with discrete spatial information (Figures 13.3c and 13.3d).
• Polyclass is outperformed by both the empirical and the fully Bayesian approach and, therefore, results are separately displayed in Figure 13.3e together with REML estimates based on two-dimensional P-splines. Presumably, the poor performance of polyclass is mainly caused by the special choice of linear splines, resulting in rather peaked estimates. Smoother basis functions, e. g. truncated cubic polynomials might improve the fit substantially but are not available in the implementation.
• Empirical and fully Bayesian estimates lead to comparable bias for the nonpara-metric effects. Results obtained with polyclass are less biased but show some spikes caused by the modeling with linear splines. Therefore, we can conclude that the poor performance of polyclass in terms of MSE is mainly introduced by additional
variance compared to empirical and fully Bayesian estimates (Figure 13.4).
• For spatial effects, both empirical and fully Bayesian estimates tend to oversmooth the data, i. e. estimates are too small for high values of the spatial functions and vice versa. In contrast, polyclass leads to estimates which are too wiggly and, therefore, overestimate spatial effects (Figures 13.5 and 13.6).
• For some simulation runs with spatial effects modeled by MRFs, no convergence of the REML algorithm could be achieved. This also happened when the spatial effect was modeled by a two-dimensional P-spline but in a much smaller number of cases. Obviously, the same convergence problems as described in Section 9 appear in a categorical setting. However, the arguments given there still hold and, hence, estimates obtained from the final (100-th) iteration are considered in these cases.
−6−5−4−3−2−1log(MSE)
MRF 2dP GRF
REML MCMC REML MCMC REML
(a) f1(x) (cat. 1)
−5−4−3−2−1log(MSE)
MRF 2dP GRF
REML MCMC REML MCMC REML
(b) f1(x) (cat. 2)
−5−4−3−2−1log(MSE)
MRF 2dP GRF
REML MCMC REML MCMC REML
(c) f2(s) (cat. 1)
−5−4−3−2−1log(MSE)
MRF 2dP GRF
REML MCMC REML MCMC REML
(d) f2(s) (cat. 2)
−8−4048log(MSE)
f1(x) (cat. 1) f1(x) (cat. 2) f2(s) (cat. 1) f2(s) (cat. 2) REML polyclass REML polyclass REML polyclass REML polyclass
(e) polyclass
Figure 13.3: Comparison of different modeling approaches: Boxplots of log(MSE).
−1−.50.51
0 .25 .5 .75 1
REML: f1(x) (cat.1)
−1−.50.51
0 .25 .5 .75 1
REML: f1(x) (cat.2)
−1−.50.51
0 .25 .5 .75 1
MCMC: f1(x) (cat.1)
−1−.50.51
0 .25 .5 .75 1
MCMC: f1(x) (cat.2)
−1−.50.51
0 .25 .5 .75 1
polyclass: f1(x) (cat.1)
−1−.50.51
0 .25 .5 .75 1
polyclass: f1(x) (cat.2)
Figure 13.4: Comparison of different modeling approaches: Bias of nonparametric esti-mates. Estimates are displayed as solid lines, the true functions are included as dashed lines.
-1.7 0 1.4
REML: f2(s) (cat. 1)
-0.6 0 0.6
REML: f2(s) (cat. 1)
-1.7 0 1.4
MCMC: f2(s) (cat. 1)
-0.6 0 0.6
MCMC: f2(s) (cat. 1)
-1.7 0 1.4
polyclass: f2(s) (cat. 1)
-0.6 0 0.6
polyclass: f2(s) (cat. 1)
Figure 13.5: Comparison of different modeling approaches: Average estimates (left panel) and empirical bias (right panel) of estimates for f2(1)(s).
-1.7 0 1.4
REML: f2(s) (cat. 2)
-0.3 0 0.3
REML: f2(s) (cat. 2)
-1.7 0 1.4
MCMC: f2(s) (cat. 2)
-0.3 0 0.3
MCMC: f2(s) (cat. 2)
-1.7 0 1.4
polyclass: f2(s) (cat. 2)
-0.3 0 0.3
polyclass: f2(s) (cat. 2)
Figure 13.6: Comparison of different modeling approaches: Average estimates (left panel) and empirical bias (right panel) of estimates for f2(2)(s).