## Abstract

Data for this research come from time series of monthly average temperatures from 28 sites over the Valle del Cauca of Colombia in South America, collected over the period 1971–2002. Because of the geographical location of the study area, monthly average temperature is affected by altitude and El Niño–La Niña (El Niño–Southern Oscillation, or ENSO phenomenon). Time series for some of the sites show a tendency to increase. Also, because of the two dry and wet periods in the study area, a seasonal pattern of behavior in monthly average temperature is seen. Linear mixed models are formulated and fitted to account for within- and between-site variations. The ENSO phenomenon is modeled by the Southern Oscillation index (SOI) and dummy variables. Spatial and temporal covariance structures in the errors are modeled individually using isotropic variogram models. The fitted models demonstrate the influence of the ENSO phenomenon on monthly average temperatures; this is seen in the maps produced from the models for ENSO and normal conditions. These maps show the predicted spatial patterns for differences in temperature throughout the study area.

## 1. Introduction

Climatic and soil factors influence crop growth and help determine the distribution of species and their adaptations to changing conditions (Atherton and Rudich 1986; Jones 2000; Villareal 1980). One of the goals of the Grupo de Investigación en Estadística Aplicada Inferir (INFERIR, the Applied Statistics Research Group at the Universidad del Valle, Cali, Colombia) is the comprehensive estimation of climatic and soil variables throughout Colombia, even at locations where these variables are not measured directly. This research is concerned with this general objective and is focused on modeling monthly average temperature in Valle del Cauca. Valle del Cauca (Fig. 1) is located in an intertropical region where monthly average temperature is affected by altitude and the ENSO phenomenon (Pabón et al. 2002).

In modeling monthly average temperature, several authors used multiple linear regression and kriging methods (Isaaks and Srivastava 1989; Cressie 1993). Spatial dependence among the data is considered in most of these studies; authors like Ayalus and Paz-González (2004) and Hoff (2001) used kriging methods, based on models for spatial covariance structures such as the spherical and the Gaussian models. Other authors like Perry and Hollis (2005a,b) used inverse distance weighting, a method that gives more weight to the closest samples and less to those farthest away from the point being estimated (Isaaks and Srivastava 1989). Lennon and Turner (1995) concluded that the best method for predicting the monthly average temperature is a thin plate spline (New et al. 1999). Pabón et al. (2002) included 400 meteorological stations in the modeling of monthly average temperature. They fitted a simple linear regression model that included altitude as a covariate. Thomas and Herzfeld (2004) obtained monthly climate maps of temperature, precipitation, and annual potential evapotranspiration, estimating values through stepwise multivariate regression, geostatistical analysis, and interpolation (the kriging method). Ninyerola et al. (2000) modeled mean temperature using regression models where the models obtained are stratified by month. Handcock and Wallis (1994) developed spatial–temporal models for average winter temperature over a region in the northern United States covering eastern Montana through the Dakotas and northern Nebraska up to the Canadian border. Other authors, like Shabbar and Barnston (1996), used canonical correlation analysis to predict the mean monthly temperature and Antonic et al. (2001) used neural networks in the modeling of seven climatic variables, including monthly mean air temperature. Andrade-Bejarano (2008) fitted a random intercept model with altitude as the covariate. The author fitted models by month and geographical position. In the research presented here, we developed and fitted linear mixed models to accommodate within- and between-site variations. We included the covariates altitude, geographical position, Southern Oscillation index, and dummy variables to model the ENSO phenomenon. We explored and modeled the spatial and temporal structures of the data by using isotropic models like the Gaussian and exponential models. The remainder of this paper is organized as follows. Section 2 describes the characteristics of the data of the case study that motivated this paper. Section 3 presents the statistical modeling of the monthly average temperature. Section 4 contains our results and discussion, and section 5 provides our conclusions and plans for future work.

## 2. Case study

The department of Valle del Cauca (study area) is located in southwest Colombia, with coordinates from 3°05′ to 5°01′N and from 75°42′ to 77°33′W. The geographical valley is limited by two mountain ranges: the Central Mountain Range and the Western Mountain Range. The Cauca River runs through Valle del Cauca. The department of Valle del Cauca has a western boundary with the Pacific Ocean. Valle del Cauca covers 2 214 000 ha, 304 379 ha (13.75%) of which is used for agriculture (Ministerio de Agricultura y Desarrollo Rural 2009). Agriculture is the main economic activity in the regional economy.

We analyzed the monthly average temperatures recorded at a network of meteorological stations (hereafter sites) in Valle del Cauca (Fig. 2) during the period 1971–2002. The sites are located in the intertropical region, from 3°19′ to 4°44′N and from 75°49′ to 76°45′W, and at altitudes of 920–1950 m above sea level. Fifteen sites are located in the valley, at altitudes up to 1100 m, with 13 in the mountains, at altitudes above 1233 m.

In the data, two sources of variation can be identified: one related to the temporal aspects of the records and one to the spatial aspects (locations of the sites).

In the aforementioned, four aspects are distinguished:

Time series of monthly average temperature over some sites show an increasing trend (Fig. 3).

The monthly average temperature shows a seasonal variation. Pabón et al. (2002) stated that in the study area there are two dry periods (January–February and July–August) and two wet periods (April–May and October–November). Figure 4 shows the seasonal variation in the study area.

The ENSO phenomenon, where the El Niño is associated with elevated monthly average temperatures, and the La Niña with reduced averages (Fig. 5).

We used time series from 28 sites and data from these can show temporal correlation, because a model for a time series will generally reflect the fact that observations close together in time will be more closely related than observations farther apart, and it is possible to find a dependency between these close observations (Box et al. 1994; Vandaele 1983).

In the spatial variation, two aspects can be identified:

Due to the altitude, monthly average temperature behaves differently in the valley and mountains. The monthly average temperature tends to be lower at higher altitudes (McGregor and Niewolt 1998; Pabón et al. 2002).

The sites within close proximity to one another are more likely to have similar values (and patterns of values) than sites farther apart.

The sources of variability motivated us to choose linear mixed models (Henderson 1982). In these models, years and sites are related with random effects and altitude, while indicators of ENSO and the geographical position are related with fixed effects.

## 3. Statistical modeling

This section describes the statistical modeling of monthly average temperature through linear mixed models, the exploration of the spatial and temporal correlations, the modeling of the spatial and temporal covariance structures in the errors, and the selection of the best of these models.

### a. The linear mixed model

In general, the linear mixed model (Laird and Ware 1982) has the form

where **y**_{i} is the *n _{i}* × 1 response vector per subject (individual or cluster)

*i*,

_{i}is an

*n*×

_{i}*p*matrix of the explanatory variables (covariates),

**is the corresponding**

*β**p*× 1 vector of regression parameters,

_{i}is an

*n*×

_{i}*q*matrix associated with random effects,

**b**

_{i}is a

*q*× 1 vector of cluster-level random effects, and

*ε*_{i}is an

*n*× 1 vector of elementary-level random errors. Matrices

_{i}_{i}and

_{i}are completely observed, and at the outset we assume that so are vectors

**y**

_{i}. The total number of observations is

*N*=

*n*

_{1}+ +

*n*, where

_{m}*m*is the number of clusters (sites). The

*N*×

*p*regression design matrix is formed by the vertical stacking of matrices

_{i},

*i*= 1, …,

*m*. Variation design matrix is composed of diagonal blocks

_{i},

*i*= 1, …,

*m*.

### b. Models for monthly average temperatures

The monthly average temperature data of this research are classified as cross-classification data (Goldstein 1995), which means that data belong simultaneously to years and sites. We fitted four linear mixed models as in (1), *M* = 1, …, 4, to the monthly average temperature data. To eliminate the seasonal effects, models were fitted for each month. Fitting models by month is useful to farmers in the study area because they make decisions on when to sow crops based on climatic information for each month. Therefore, we studied 4 × 12 = 48 model fits.

The model corresponding to the *j*th month, *j* = 1, …, 12, is

where *i* = 1, …, 27 sites and *k* = 1, …, 32 yr.

All four models have the same random part. The difference among models is in the fixed part. The set of variables for *M _{h}*,

*h*= 1, …, 4 for the fixed part of each model includes

where the index 0 corresponds to the intercept (*x*_{0ik} ≡ 1). The variables 1–9 and its parameters *β* are listed in Table 1.

#### 1) Random part of the model in (2)

a random effect for years (

*b*_{2k}),*k*= 1, …, 32 yr, which models the year-to-year variation in monthly average temperature due to the ENSO phenomenon, as well as other factors;random coefficients to model the trends of the time series, as described in section 2, include the site (intercept) (

*b*_{1i}) and site year effect (slope) (*b*_{3i}), which are used to model separate regressions of monthly average temperature over time for each site (*i*= 1, …, 27 sites); note that because straight lines are fitted, the site effects (intercepts) and site year effects (slopes) are correlated within each site;a vector of the random errors

*ɛ*; and_{ik}as indicated before, the variation design matrix is composed of diagonal blocks

_{i},*i*= 1, …,*m.*

#### 2) Fixed part of the model in (2)

As indicated in section 2, the valley and mountains behave differently with altitude (m). Furthermore, average temperature decreases by 0.625°C for each 100-m increase in altitude (Pabón et al. 2002) and these differences are particularly significant in the tropics (McGregor and Niewolt 1998).

As was outlined in sections 1 and 2, the El Niño and La Niña phenomena affected the monthly average temperature values, and the Southern Oscillation index (SOI) was included in the model in (2) as a measurement of these phenomena.

The SOI with lag two (SOI_LAG2) from the previous two months was also included. Madl (2000) indicated the influence of the months previous to the El Niño phenomenon on the weather, as the author stated that in the months preceding an El Niño event, the normal weather pattern breaks down. For some reasons not yet well understood, the westward atmospheric pressure gradient decreases.

Due to the monthly average temperature behaving differently in the valley and the mountains, the geographical position also was included as a dummy variable in the models. The monthly average temperature oscillates between 15.2° and 23.8°C in the mountains, and between 21.6° and 29.9°C in the valley.

We also modeled the El Niño and La Niña phenomena for the current month and for the previous two lag months through dummy variables as follows:

for the current month's El Niño,

*x*_{5k}= 1,*x*_{6k}= 0; La Niña,*x*_{5k}= 0,*x*_{6k}= 0; and for normal conditions,*x*_{5k}= 0,*x*_{6k}= 1; andfor El Niño in the two lag months,

*x*_{7k}= 1,*x*_{8k}= 0; La Niña in the two lag months,*x*_{7k}= 0,*x*_{8k}= 0; and normal conditions in the two lag months,*x*_{7k}= 0,*x*_{8k}= 1.We included the year translated (year minus 2002) in the random part of the models to model the random year site slopes. This variable should also be included in the fixed part of the models, because of the statistical specification of the random coefficients (Verbeke and Molenberghs 2000).

The random effects *b*_{1i} and *b*_{3i} are correlated and have the following distribution:

, and *ɛ _{ik}* ~

*N*(0,

*σ*

^{2});

*b*

_{2k}and

*ɛ*are mutually independent random samples.

_{ik}The spatial and temporal correlated error models are given by different *M _{h}*, for

*h*= 1, …, 4. For the vector of the

*ik*random term

*ɛ*, we assumed that

_{ik}*ɛ*~

*N*(0,

*σ*

^{2}

*F*), where

*F*is an isotropic covariance function (Chilès and Delfiner 1999). These functions are specified in section 3d.

We fitted the four models by month using the Proc Mixed of SAS version 9.1 (SAS 2004) and the restricted maximum likelihood method (REML; Patterson and Thompson 1971). We tested the statistical significance of the parameter estimates in the models, using the *t* test (Verbeke and Molenberghs 2000). The covariance structure used for year are variance components, and unstructured is used for the random intercept and slope. Initially, we fitted the models assuming independent errors. Some observations from one site (Zaragoza) showed large residuals so we developed a simulation study to decide whether these observations were outliers (Andrade-Bejarano and Longford 2010). The results of the simulation study showed that these observations were indeed outliers, thus eliminated this site from our study. The most likely explanation for the inconsistent data from Zaragoza is failure in the measurement equipment. In a second step, we explored the spatial and temporal correlations in the errors.

### c. Exploration of the spatial and temporal correlation in the errors

We explored the spatial and temporal correlations of the errors by using variograms. We also used autocorrelation function graphs of time series from sites, to evaluate the temporal correlation in the errors. The variogram of residuals for the spatial correlation in the errors is

where *N*(*h*) is the count of the pairs of points (sites) separated (approximately) by the lag *h* and *r*(*s _{ik}*) and are the residual values of the variable in the point

*s*and , respectively, for year

_{i}*k*. In the case of the temporal variograms,

*r*(

*s*) and are replaced by

_{ik}*r*(

*t*) and , where

_{ik}*r*(

*t*) corresponds to the residual value in time

_{ik}*k*and to the residual value in time

*k*′ for site

*i*.

Finally, we worked with 27 sites (as was indicated before, one of them was removed from the analysis as a result of the outlier analysis). These sites are located in the valley and in the mountains. The data span a 32-yr period; however, in some sites and years data are missing. Due to the missing data, we calculated a variogram weighted by the number of pairs of points we used to calculate each variogram per year and, in the case of the temporal variograms, the number of pairs of points per site, yielding

where *γ*_{1}(*h*), …, *γ _{n}*(

*h*) are the variograms from

*n*years and

*N*

_{1}(

*h*), …,

*N*(

_{n}*h*) is the count of pairs of points separated (approximately) by the lag

*h*from each year, or sites in the case of temporal variograms.

A common form for choosing the lag *h* is to define a set of distance intervals and to consider all the pairs of points whose distances are in the same interval (Samper and Carrera 1990). We followed this criterion and in the case of the spatial variograms we chose a separation length *h* (lag distance) of 10 km. In the case of temporal variograms, we used a lag *h* of 1 yr.

Because of the small number of sites (27 in total), we did not explore anisotropy. Authors such as Hoff (2001) also did not explore anisotropy, despite the fact that he used 190 sites for the period 1979–93.

### d. Modeling of the spatial and temporal correlations in the errors and choosing the best model

After we explored the spatial and temporal correlations in the errors by variograms, we fitted the linear mixed models and modeled the spatial and temporal covariance structures in the errors by using the three isotropic variogram models: Gaussian, exponential, and spherical (Cressie 1993). In the modeling of the temporal covariance structure of the errors, we also used the autoregressive model [AR(1); Verbeke and Molenberghs (2000)]. For the spatial correlation in the errors, Gaussian, exponential, and spherical models are defined by

Gaussian model,

exponential model,

- spherical model, where is the distance between sites
*i*and*i*′ and*ρ*_{s}is the range of the spatial variogram. In the case of the temporal correlation of the errors, is replaced by , meaning the lag time and*ρ*is replaced by_{s}*ρ*, the range of the temporal variogram (Chilès and Delfiner 1999)._{t}

We determine the best model through the likelihood ratio test (Verbeke and Molenberghs 2000), which compares models with independent errors to models with spatial or temporal correlated errors. The likelihood ratio test (LRT) can be expressed as (Verbeke and Molenberghs 2000):

where *l*_{REML0} is the logarithm from the maximum likelihood estimation computed under the null model using the REML method and *l*_{REML1} is the logarithm from the maximum likelihood estimations computed under the alternative model. In the case of this research, the null model is the model with independent errors and the alternative model is the model with spatial or temporal covariance structures in the errors, as modeled by the spherical, Gaussian, or exponential models specified in section 3d. The model with the largest LRT value was selected as the best model.

For a given set of fixed effects in the models with covariates *M _{h}*,

*h*= 1, …, 4, we tested the hypotheses of no spatial and no temporal correlations. The null hypotheses are

*ρ*

_{s}= 0 and

*ρ*

_{t}= 0. To determine the statistical significance of these hypotheses, we derived by simulation the null distribution of LRT and compare the calculated values in (5) with the 95th percentile of this derived null distribution of the LRT (Andrade-Bejarano 2009).

## 4. Results and discussion

### a. Results of the exploration of the spatial and temporal correlations of the errors

For all models except August, the omnidirectional spatial variograms by month generally show a spatial correlation up to 21.2 km [Fig. 6 illustrates it, for the first model specified in Eq. (2), Table 2, section 3b]. Beyond this distance no months or models evaluated showed any remaining correlation. The variograms for August did not show a spatial correlation: the initial value of *γ* is 0.27, which declined at 10.61 km and at 30.6 km *γ* increases again.

Isaaks and Srivastava (1989) defined a nugget effect as “the vertical jump to the value of 0 at the origin to the value of the variogram at extremely small separation distances.” However, despite the fact that the value of the variogram for *h* = 0 is strictly 0, several factors, such as sampling error and short-scale variability, may cause sample values separated by extremely small distances to be quite dissimilar, thus generating a nugget effect. In spatial variograms for monthly average temperature, we did not find a nugget effect, as the distance between the closest meteorological stations is 1.84 km and they showed similar monthly average temperature values.

For all months and fitted models, variograms show a temporal correlation in the errors (Fig. 7 shows this for the first model). Temporal variograms also do not show a nugget effect. The lag distance is 1 yr, and thus a possible explanation is that a 1-yr time lag may be too small to find a nugget effect.

As we indicated above, due to the small number of sites we could not explore anisotropy. We assume isotropy (Chilès and Delfiner 1999) and we modeled the covariance structure of the errors using isotropic covariance models. We eliminated the seasonal component of the time series by fitting models by month. The trends found in some time series were modeled through a random intercept and the random slope, which are included in the four fitted models.

We separately analyzed the spatial and temporal correlations in the errors. A spatiotemporal covariance structure in the errors was not developed. An initial exploration of the spatiotemporal variograms of monthly average temperature was carried out, following the criteria of Cressie and Huang (1999). Spatial variograms for the first time lag show a concave form and spatial continuity behavior (Andrade-Bejarano 2009).

Authors like Cressie (1993), Gneiting (2002), Ma (2003a,b), Kolovos et al. (2004), and Stein (2005) have made contributions to the derivation of spatiotemporal covariance functions. Banerjee et al. (2004) discuss the nonseparable stationary covariance functions developed by Cressie (1993), Gneiting (2002), and Stein (2005); they indicate the difficulties in the computation of cov(*s* − *s*′, *t* − *t*′). Similarly, Stein (2005) states that it is necessary to develop theoretical frameworks for the differences in the models and the methods for model fitting in order to compare models and diagnostics and to assess the adequacy of the models.

### b. Choosing the best model

As mentioned in section 3d, we determined the best model to model the spatial and temporal covariance structures in the errors using the LRT given in Eq. (5). Table 2 shows the LRT values from the first model for the mathematical difference between twice the *l*_{REML0} value of models with independent errors (1) and *l*_{REML1} values of models with spatial or temporal covariance structures in the errors, modeled by (a) the spherical model in (2), (b) the Gaussian model in (3), and (c) exponential model in (4).

The Gaussian model is the best scheme for modeling the spatial covariance structure in the errors, because in most of the months it shows the largest LRT values among the three models evaluated. This result is outlined in Table 2 as LRT(1–3). Figure 6 shows the proximity of the omnidirectional variograms to the Gaussian model. For August, the Newton–Rhapson algorithm for the estimation of LRT_{1} values for the spherical model did not converge (Table 2). For this month, we also obtained zero LRT values. These results are due to there being no spatial correlation in the errors for this month for any models, as seen in the variogram (Fig. 6).

In the case of the temporal covariance structure in the errors, the exponential model and the AR(1) model display the largest values of likelihood ratio test and also the same values (Andrade-Bejarano 2009). Verbeke and Molenberghs (2000) and Avalos et al. (2008) also found the Exponential model to be a good model for the temporal covariance structure in the errors. Figure 7 shows the proximity of the omnidirectional variograms to the exponential model. As a consequence, we chose the exponential model to model the temporal covariance structure of the errors.

### c. Results of the random part estimates of the models

For all months and fitted models (assuming independent errors or spatial and temporal covariance structures in the errors), variance related to site shows the largest values among the covariance parameter estimates, and exhibits the lowest values. Table 3 displays the results for months during the driest and the wettest periods of the study area, for the first model. For the first, second, and fourth models [section 3b, Eq. (2), Table 2], the largest estimated value of the variance due to years , and the standard error and interval values, correspond to January (Table 3), and for the third model to February (Andrade-Bejarano 2009). We obtained similar results in models assuming independent errors. January and February are the driest months in the study zone and temperature increases during this period, especially when the El Niño phenomenon occurs (Pabón et al. 2002).

We found a high correlation between *b*_{1} and *b*_{3}, where the values fluctuate between 0.69 and 0.95. This result is expected, due to the assumption of correlation between the random intercept and the random slope. For all the models, July has the largest estimated value, standard error, and interval values for , while October has the lowest estimated values (Table 3). We obtained the same result in models where we assumed independent errors.

For all the models with spatial covariance structure in the errors, estimated *ρ* values, which are the ranges of the variograms, fluctuate between 2.6 and 5.3 (Andrade-Bejarano 2009), with May and December being the months with the lowest estimated values and the highest standard error values. For the models with a temporal covariance structure in the errors, estimated *ρ* values (the ranges of the variograms) fluctuate between 0.7 and 1.6.

### d. Results of the fixed part estimates of the models

#### 1) Intercept

Estimated values of *β*_{0}, the intercept parameter, are significant at *α* = 0.05 for all months and models with spatial and temporal covariance structures in the errors. Tables 4 and 5 display the results for the months for the driest and the wettest periods of the study area for the four models. We obtained similar results with models assuming independent errors.

#### 2) Altitude variable

For all months and models assuming independent errors and with spatial and temporal covariance structures in the errors, the estimated values of *β*_{1} are statistically significant at *α* = 0.05 (Tables 4 and 5), which indicate that the altitude variable contributes to the explanation of the monthly average temperature. This result agrees with Pabón et al. (2002) regarding the influence of altitude on temperature. Florio et al. (2004), Perry and Hollis (2005a,b), Thomas and Herzfeld (2004), Ayalos and Paz-González (2004), Pabón et al. (2002), Hoff (2001), Ninyerola et al. (2000), Handcock and Wallis (1994), and Antonic et al. (2001) also included altitude as a covariate in temperature modeling.

#### 3) Geographical position variable

In most months and models assuming independent errors and with spatial and temporal covariance structures in the errors, the estimated *β*_{9} values are statistically significant at *α* = 0.05 (Tables 4 and 5). This indicates that the geographical position contributes to the explanation of the monthly average temperature. However, for all fitted models, estimated values of *β*_{9} for May, October, and November are not statistically significant. These months belong to the wet period during the study area.

#### 4) SOI and SOI_LAG2 variables

As indicated in section 3b, SOI and SOI_LAG2 were included in the first and second models, respectively. For the first model, assuming independent errors and with spatial and temporal covariance structures in the errors, the estimated *β*_{2} values are statistically significant at *α* = 0.05 (Tables 4 and 5), except for May.

In the second model, assuming independent errors and with spatial and temporal covariance structures in the errors, the estimated *β*_{4} values are statistically significant at *α* = 0.05 (Tables 4 and 5). For all months, these results indicate that the Southern Oscillation index for the current month and the Southern Oscillation index of 2 months prior explain the monthly average temperature in the study area.

#### 5) Dummy variables included in the modeling of the El Niño and La Niña phenomena

As indicated in section 3b, we used two dummy variables, with parameters *β*_{5} and *β*_{6}, for the modeling of the El Niño and La Niña phenomena for the current month, in the third model. Similarly, we used two other dummy variables for the modeling of the El Niño and La Niña phenomena 2 months prior, with parameters *β*_{7} and *β*_{8}, in the fourth model.

In the third model, assuming independent errors and with spatial and temporal covariance structures in the errors, the estimated *β*_{5} parameters are statistically significant at *α* = 0.05 (Tables 4 and 5) during most months. The estimated *β*_{6} parameters are statistically significant at *α* = 0.05 for 6 months.

In the fourth model assuming independent errors and with spatial and temporal covariance structures in the errors, the estimated *β*_{7} parameter are statistically significant at *α* = 0.05. However, *β*_{8} are not statistically significant in most months.

#### 6) Year translated

For all models with independent errors, and with spatial and temporal covariance structures in the errors, *β*_{3}, the parameter of the year translated, is statistically significant at *α* = 0.05 for October, November (Tables 4 and 5), and December (Andrade-Bejarano 2009). For May and July, in most models with temporal covariance structure in the errors, the estimated *β*_{3} parameters are statistically significant at *α* = 0.05.

#### 7) Statistical significance of the spatial and temporal correlation in the errors

The results show that for the four models, except for May, August, and December, there is a statistically significant spatial correlation in the errors at *α* = 0.05 (Andrade-Bejarano 2009). For all the models and months there is a statistically significant temporal correlation in the errors at *α* = 0.05. The likelihood ratio test values computed by comparing models with independent errors and models with temporal covariance structures in the errors exhibit high values. Because we found statistically significant spatial and temporal correlations in the errors, we needed to model them, in this case, through the Gaussian and exponential models, for the spatial and temporal covariances structures of the errors, respectively.

### e. Predictions using fitted models

The differences among the models were described in section 3b. All the models have the same random part. The difference between the first and second models is in the inclusion of SOI and SOI_LAG2, respectively. The difference between the third and fourth models is the inclusion of the El Niño and la Niña phenomena for the current month and two lag months, respectively, using dummy variables in both cases. Because of the statistical similarities in the fixed part between the first and second models and between the third and fourth models, respectively, we used the first and the third models, with spatial and temporal covariance structures in the errors, to obtain the predictions.

Maps of predicted values obtained from the first model (Fig. 8) and third model (Andrade-Bejarano (2009), Fig. 11.3) show similar patterns: during the La Niña phenomenon [Fig. 8a and Andrade-Bejarano (2009, Fig. 11.3(a)], predicted monthly average temperature values vary between 24° and 16.5°C and the maps show a decrease in these values with respect to normal conditions, with values between 24.5° and 17.2°C [Fig. 8c and Andrade-Bejarano (2009), Fig. 11.3(c)]. Maps for the two models also show that during the El Niño phenomenon [Fig. 8b and [Andrade-Bejarano (2009), Fig. 11.3(b)], predicted monthly average temperature values increase with respect to normal conditions with values between 25.3° and 18°C.

In February, the first model gives, on average, a difference in the mean temperatures for El Niño conditions of 0.9°C above normal, whereas La Niña shows a decrease of 0.5°C. The third model gives, on average, a difference in the mean temperatures for El Niño conditions of 0.7°C above normal, whereas La Niña shows a decrease of 0.4°C.

Maps of predicted values of monthly average temperature obtained from the first and third models with temporal covariance structure in the errors (Andrade-Bejarano 2009, Figs. 11.2 and 11.4), for February, show very similar patterns to maps obtained by the models with spatial covariance structure in the errors. Under normal conditions the values vary between 24.5° and 17.2°C. During the La Niña phenomenon, the predicted monthly average temperature decreases, and during the El Niño it increases. These results agree with Pabón et al. (2002) concerning the increase and decrease of the mean temperature values during the El Niño and La Niña phenomena, respectively.

The differences in mean temperatures, between El Niño conditions and normal and between La Niña conditions and normal, are equal to the differences obtained in the first and third models with spatial covariance structures in the errors.

We can see a very similar pattern between maps from the first (Fig. 9) and third models (Andrade-Bejarano 2009, Fig. 11.7) with spatial covariance structures in the errors, for November. These maps show a decrease in mean monthly average temperature values during the La Niña phenomenon (values vary between 23.5° and 16.5°C) in comparison to normal conditions (varying between 23.8° and 17.0°C). Maps also show an increase in the mean monthly average temperature values (values vary between 24° and 17.1°C) during the El Niño phenomenon in comparison to normal conditions. The first model gives, on average, a difference in the mean temperatures for El Niño conditions of 0.2°C above normal, whereas La Niña shows a decrease of 0.1°C. The third model gives, on average, a difference in the mean temperatures for El Niño conditions of 0.2°C above normal, whereas La Niña shows a decrease of 0.3°C.

Maps for models with temporal covariance structures in the errors for November (Andrade-Bejarano 2009, Figs. 11.6 and 11.8) show a very similar pattern to maps obtained by models with spatial covariance structures in the error. The differences in mean temperatures, between El Niño conditions and normal and between La Niña conditions and normal, are equal to the differences obtained in the first and third models with spatial covariance structures in the errors.

In February, the first model predicts larger temperature increases (or decreases) of El Niño and La Niña from normal: the largest changes are El Niño increases of 0.9° and 0.7°C above normal for the first and third models. La Niña shows decreases of 0.5° and 0.4°C from normal, for the two models, respectively.

For November, in the first model El Niño is 0.2°C above normal and La Niña is 0.1°C below normal; meanwhile, for the third model El Niño is 0.2°C above normal and La Niña is 0.3°C below normal.

These results show an important difference in the effects of the El Niño and La Niña phenomena in February and November. The effect increases temperatures by 0.9° and 0.7°C in February, for the two models, and but only 0.2°C in November. However, the effect of La Niña is relatively uniform, with decreases of 0.5° and 0.4°C in February for the first and third model, respectively, and 0.1° and 0.3°C in November, for the first and third models.

The average patterns of change between El Niño, normal, and La Niña conditions, as shown on these maps, are consistent with respect to the first and third models with spatial and temporal covariance structures in the errors. The first model uses the quantitative value of the SOI, which is a measure of the strength of the El Niño and La Niña phenomena, whereas the third model is based on the classification of years and does not include a measure of strength.

## 5. Conclusions

We fitted four models through the linear mixed models framework. We included the same random effects in our four models but each with different fixed parts. The trends of some time series are well modeled by the inclusion of the site effect, which provides the intercept term for each site, and a site time interaction, which provides the slopes for each site in the models. The random effect due to year shows variability in the estimated values among months; the values for the driest period (January–February) are larger than values for the wettest period (October–November). This is reasonable given the increase in temperature during the driest period (Pabón et al. 2002) and the decrease in temperature during the wettest period, as well as the substantial influence of the El Niño and the La Niña phenomena during the 1971–2002 period.

All the models show the influence of the altitude variable upon the monthly average temperature, which is in agreement with the physical phenomenon due to altitude. The geographical position in the valley and the mountains also explains the monthly average temperature in all the models and during most months.

We modeled the El Niño and La Niña phenomena by using the Southern Oscillation index (SOI) for the current month and the Southern Oscillation index for the two previous months (SOI_LAG2). Both models that include these variables (the first and second models) show the influence of the El Niño and La Niña phenomena in modeling the monthly average temperature. We also modeled these phenomena by using dummy variables and we found equivalent results among the maps of predicted values obtained with the first and third models, in the sense that both models properly model the climatic conditions in the study area, and they capture the influence of the El Niño and La Niña phenomena on the monthly average temperature, showing increases during El Niño and decreases during La Niña. The maps also reflect the behavior of the monthly average temperature under normal conditions and the increases and decreases of the monthly average temperature during the months belonging to the driest and the wettest periods, respectively. In February, the first model predicts larger than normal temperature increases (or decreases) in El Niño and La Niña compared with the third model.

We separately analyzed the spatial and the temporal correlations in the errors. A spatiotemporal covariance structure in the errors was not developed. An initial exploration of the spatiotemporal variograms of monthly average temperature was carried out. Spatial variograms for the first time lag show a concave form and spatial continuity behavior.

As was indicated previously, we used linear mixed models and this framework has not been used to model the monthly average temperature. Other research into the modeling of monthly average temperature has been conducted within the following frameworks: 1) regression and interpolation methods and 2) canonical analysis and neural networks. The inclusion of the random effects (sources of variability), which in our case are the meteorological station and year, provides better estimated *β* values (i.e., correct signs and reduction of the standard errors) in comparison with estimated *β* values obtained by ordinary least squares method (Demidenko 2004; Verbeke and Molenberghs 2000; Brown and Prescott 1999). We also included covariables to model the influence of the El Niño and La Niña phenomena. In the literature, this influence of the ENSO phenomenon has not been included in the models fitted by other authors.

Future work is recommended to develop the spatiotemporal covariance structure in the errors for models fitted in this research. This promises to be a challenge in spatiotemporal modeling; however, as Stein (2005) has indicated, many theoretical aspects need to be researched, including ways to estimate the spatiotemporal covariance structure. Another suggestion for future work is the modeling of temperature for all of Colombia. Colombian geography is very complex, as three mountains ranges split from the Andes across the country. Similarly, the country has shores on two oceans (Pacific and Atlantic). Furthermore, Colombian ecosystems consist of savanna, desert, jungle, and rain forest, generating a diversity of climates. In accordance with the results of the present research, linear mixed models, geostatistics, and time series analysis are promising statistical frameworks for such a project.

## Acknowledgments

This work was sponsored by Universidad del Valle, Cali, Colombia. Data for this research were supplied by Instituto de Hidrología, Meteorología y Estudios Ambientales de Colombia (IDEAM) and Corporación Autónoma Regional del Valle del Cauca, Colombia (CVC).

## REFERENCES

*GeoENV VI—Geostatistics for Environmental Applications,*A. Soares, M. Pereira, and R. Dimitrakopoulos, Eds., Quantitative Geology and Geostatistics, Vol. 15, Springer Science+Business Media B.V., 361–371.

*The Tomato Crop: A Scientific Basis for Improvement*. Chapman and Hall, 684 pp.

*N*

_{2}

*O*fluxes from differently treated soils.

*GeoENV VI—Geostatistics for Environmental Applications,*A. Soares, M. Pereira, and R. Dimitrakopoulos, Eds., Quantitative Geology and Geostatistics, Vol. 15, Springer Science+Business Media B.V., 361–371.

*GeoENV IV: Geostatistics for Environmental Applications,*X. Sanchez-Vila, J. Carrera, and J. Gómez-Hernández, Eds., Kluwer Academic, 520–521.

*Hierarchical Modelling and Analysis for Spatial Data*. Chapman and Hall, 452 pp.

*Time Series Analysis, Forecasting and Control*. Prentice- Hall, 598 pp.

*Applied Mixed Models in Medicine*. John Wiley and Sons, 408 pp.

*Geostatistics: Modeling Spatial Uncertainty*. John Wiley and Sons, 695 pp.

*Statistics for Spatial Data*. John Wiley and Sons, 900 pp.

*Mixed Models: Theory and Applications*. Wiley Series in Probability and Statistics, John Wiley and Sons, 704 pp.

*Multilevel Statistical Models*. 2nd ed. Kendall's Library of Statistics, Vol. 3, Oxford University Press, 178 pp.

*Valle del Cauca: Aspectos Geográficos*. La Subdirección, 148 pp.

*An Introduction to Applied Geostatistics*. Oxford University Press, 561 pp.

*Tomato Plant Culture: In the Field, Greenhouse and Home Garden*. CRC Press, 399 pp.

*Tropical Climatology: An Introduction to the Climates of the Low Altitudes*. John Wiley and Sons, 339 pp.

*La Atmósfera, el Tiempo y el Clima*. Instituto de Hidrología, Meteorología y Estudios Ambientales, Bogota, Colombia, 91 pp.

*Geoestadśtica: Aplicaciones a la Hidrogeología Subterránea*. Centro Internacional de Métodos Numéricos en Ingeniería, 484 pp.

*Applied Time Series and Box-Jenkins Models*. Academic Press, 568 pp.

*Linear Mixed Models for Longitudinal Data*. Springer Science+Business Media B.V., 568 pp.

*Tomatoes in the Tropics*. Westview Press, 174 pp.