Advanced climwin

Adding interactions with climate

Until now, we have assumed that the effects of climate do not interact with any other variables. This may not always be the case. As of v1.1.0, users can build a custom model structure with the ability to include interactions between climate and other variables. In our example using the Mass dataset, we might expect the effect of temperature on chick mass to vary with the average age of the mothers.

Firstly, it is necessary to create a variable called ‘climate’ in the biological dataset. The variable ‘climate’ doesn’t need to contain any information, it is simply a placeholder to specify where climate data will be included in the model.

NOTE: It is necessary to call the new object ‘climate’ so that climwin can distinguish this variable from others. Other variable names will not work.

Mass$climate <- 1

Now that we have our ‘climate’ variable, we can add it into our baseline model. We can then specify how we want our climatic data to interact with other variables in the model.

baseline = lm(Mass ~ climate*Age, data = Mass)

In this example, climatic data from each tested climate window will be included in the model interacting with parent age.

Interaction <- slidingwin(xvar = list(Temp = MassClimate$Temp),
                          cdate = MassClimate$Date,
                          bdate = Mass$Date,
                          baseline = lm(Mass ~ climate*Age, data = Mass),
                          cinterval = "day",
                          range = c(150, 0),
                          type = "absolute", refday = c(20, 05),
                          stat = "mean",
                          func = "lin")

If we look at the best model from this analysis, we can see the interaction between temperature and age included in the model coefficients.

summary(Interaction[[1]]$BestModel)

Call:
lm(formula = yvar ~ climate + Age + climate:Age, data = modeldat)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6266 -1.5716  0.2878  1.6086  4.7510 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 170.2628     7.1678  23.754  < 2e-16 ***
climate      -5.5466     0.9200  -6.029 3.32e-07 ***
Age          -2.6046     2.6603  -0.979    0.333    
climate:Age   0.4024     0.3395   1.185    0.242    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.449 on 43 degrees of freedom
Multiple R-squared:  0.7778,    Adjusted R-squared:  0.7623 
F-statistic: 50.17 on 3 and 43 DF,  p-value: 4.267e-14

Mixed effects models

It will often be important to fit models with random effects terms (e.g. Individual ID). climwin includes functionality for mixed effects models using both lme4 and nlme. Random effects terms can easily be added to the baseline model. We can use the ‘Offspring’ dataset in R as an example.

baseline = lmer(Offspring ~ 1 + (1|BirdID), data = Offspring, REML = F)

NOTE: Mixed effects models should be fitted using a maximum likelihood approach (i.e. REML = FALSE). If this is not specified in the baseline model, climwin* will make this adjustment internally.*

In a more specific scenario, we may want to include our climatic variable as a random effect in our model. As above, this can be done by creating a ‘climate’ variable in our dataset and including it in the baseline model.

baseline = lmer(Mass ~ 1 + (1|climate), data = Mass, REML = F)

A worked example

Below we will provide a worked example using the Mass and MassClimate dataset. In this example, lets imagine we are interested in testing the impact of the number of days above freezing on our mass response variable. To do this we would set both our ‘upper’ and ‘binary’ parameters.

upper = 0 binary = TRUE

As we are interested in measuring the number of days above freezing, we set our stat parameter to ‘sum’. Otherwise, we used model parameter values identical to our earlier vignette.

library(climwin)

MassWin <- slidingwin(xvar = list(Temp = MassClimate$Temp),
                      cdate = MassClimate$Date,
                      bdate = Mass$Date,
                      baseline = lm(Mass ~ 1, data = Mass),
                      cinterval = "day",
                      range = c(150, 0),
                      upper = 0, binary = TRUE,
                      type = "absolute", refday = c(20, 05),
                      stat = "sum",
                      func = "lin")

When we examine the best model data, we can see that our climate data is now count data.

head(MassWin[[1]]$BestModelData)

A worked example

Spatially replicated data can be analysed using the slidingwin function with the addition of the ‘spatial’ parameter. As with regular slidingwin analysis, analysis with spatial replication requires a separate biological and climate dataset. However, these datasets should now contain an additional variable which specifies the site at which biological and climate data were collected. Below, we have called this parameter ‘SiteID’.

NOTE: The climate dataset for spatially replicated climwin analysis will often include duplicate dates. In a regular climwin analysis this will lead to errors.

With these new datasets, we can carry out a slidingwin analysis with the addition of a ‘spatial’ parameter.

MassWin <- slidingwin(xvar = list(Temp = Climate$Temp),
                      cdate = Climate$Date,
                      bdate = Biol$Date,
                      baseline = lm(Mass ~ 1, data = Biol),
                      cinterval = "day",
                      range = c(150, 0),
                      type = "absolute", refday = c(20, 05),
                      stat = "mean",
                      func = "lin", spatial = list(Biol$SiteID, Climate$SiteID))

The ‘spatial’ parameter is a list item that includes the SiteID variable for the biological and climate datasets respectively. When slidingwin fits individual climate windows, each biological record will be matched with the climate data from the corresponding site.

A worked example

weightwin can often be useful to use with climate data where we have already identified a climate window using slidingwin. Here, we will use weightwin to further investigate the Mass and MassClimate data included with the climwin package.

The basic parameters in weightwin are the same as slidingwin (though note the absence of the ‘stat’ parameter). In addition however, we must designate which weight distribution we want to use. In this case we consider a Weibull function.

weightfunc = "W"

Next, we must set the location, scale and shape values for the starting distribution that will be used to begin the optimisation procedure. The default values (3, 0.2, 0) are often appropriate for fitting Weibull distributions. However, you can explore different parameter values using the explore function.

The optimisation function used in weightwin has the potential to get ‘stuck’ during the optimisation routine. With the argument ‘n’ we can run our weightwin analysis multiple times to ensure we are finding the true optimal weight distribution.

n = 5

set.seed(100)

weight <- weightwin(n = 5, xvar = list(Temp = MassClimate$Temp), cdate = MassClimate$Date, 
                    bdate = Mass$Date, 
                    baseline = lm(Mass ~ 1, data = Mass), 
                    range = c(150, 0), 
                    func = "lin", type = "absolute", 
                    refday = c(20, 5), 
                    weightfunc = "W", cinterval = "day",
                    par = c(3, 0.2, 0))

If we look at the iterations object in the weightwin output, we can see the outcome of all 5 iterations.

weight$iterations

start_shape, start_scale and start_location show the starting parameters used for each iteration of weightwin. In the deltaAICc column we can see the outcomes of the different iterations. While 4 of the iterations reached the same peak, the 5th got stuck during the optimisation routine. This is a useful demonstration of why multiple iterations are important.

It seems that there is a clear peak that our optimisation routine has been able to identify, next we can visualise that peak using the explore function.

We can use the shape, scale and location parameters of the first iteration to fit our plot.

weight[[1]]$WeightedOutput

In the above plot, we can see that the importance of temperature declines rapidly either side of the peak.

If we compare this value to the delta AICc obtained from the slidingwin function we can see that the weightwin function is better able to explain variation in our mass parameter as we have removed the assumption of a uniform weight distribution.

Pros and cons

weightwin can provide greater detail on the relationship between climate and the biological response, such as the occurrence of exponential functions. Additionally, by using more diverse weight distributions, weightwin will often generate models with better ΔAICc values, which may be especially important when users are most interested in achieving high explanatory power. Furthermore, by using an optimisation routine weightwin often tests far fewer models than slidingwin, allowing for more rapid analysis.

Despite these benefits, weightwin will not always be the most appropriate function for all scenarios. Firstly, the nature of the fitted weight distributions means that weightwin can only detect single climate signals, which forces users to detect and compare potential climate signals with separate analyses. This is not an issue using slidingwin.

Secondly, weightwin can be a more technically complex process. While the above example works easily, in other scenarios the optimisation procedures may often get stuck on false optima or fail to converge. In these cases, users may be required to test different starting parameters and adjust optimisation characteristics such as step size. Often this procedure can be inhibitive for users with less technical knowledge.

Finally, weightwin can only be used for testing mean climate, with no capacity to consider other aggregate statistics. Therefore, whether one chooses to use weightwin or slidingwin will largely depend on the summary statistic used, the level of detail desired, and the user’s technical knowledge.

To ask additional questions or report bugs please e-mail:

[email protected]

or visit our google group