IbkTemperature: Temperature Observations and Numerical Weather Predictions...
In lmSubsets: Exact linear subset regression

Description Usage Format Source References Examples

00UTC temperature observations and corresponding 24-hour reforecast ensemble means from the Global Ensemble Forecast System (GEFS, Hamill et al. 2013) for SYNOP station Innsbruck Airport (11120; 47.260, 11.357) from 2011-01-01 to 2015-12-31.

1	data("IbkTemperature")

A data frame containing 1824 daily observations/forecasts for 42 variables. The first column (temp) contains temperature observations at 00UTC (coordinated universal time) in degree Celsius, columns 2–37 are 24-hour lead time GEFS reforecast ensemble means for different variables (see below). Columns 38–42 are deterministic time trend/season patterns.

temp: Observed temperature at Innsbruck Airport.
tp: Total accumulated precipitation.
t2m: Temperature at 2 meters.
u10m: U-component of wind at 10 meters.
v10m: V-component of wind at 10 meters.
u80m: U-component of wind at 80 meters.
v80m: U-component of wind at 80 meters.
cape: Convective available potential energy.
ci: Convective inhibition.
sdlwrf: Surface downward long-wave radiation flux.
sdswrf: Surface downward short-wave radiation flux.
sulwrf: Surface upward long-wave radiation flux.
suswrf: Surface upward short-wave radiation flux.
ghf: Ground heat flux.
slhnf: Surface latent heat net flux.
sshnf: Surface sensible heat net flux.
mslp: Mean sea level pressure.
psfc: Surface pressure.
pw: Precipitable water.
vsmc: Volumetric soil moisture content.
sh2m: Specific humidity at 2 meters.
tcc: Total cloud cover.
tcic: Total column-integrated condensate.
tsfc: Skin temperature.
tmax2m: Maximum temperature.
tmin2m: Minimum temperature.
st: Soil temperature (0-10 cm below surface).
ulwrf: Upward long-wave radiation flux.
wr: Water runoff.
we: Water equivalent of accumulated snow depth.
wp: Wind mixing energy.
w850: Vertical velocity at 850 hPa surface.
t2pvu: Temperature on 2 PVU surface.
p2pvu: Pressure on 2 PVU surface.
u2pvu: U-component of wind on 2 PVU surface.
v2pvu: U-component of wind on 2 PVU surface.
pv: Potential vorticity on 320 K isentrope.
time: Time in years.
sin, cos: Sine and cosine component of annual harmonic pattern.
sin2, cos2: Sine and cosine component of bi-annual harmonic pattern.

Observations: http://www.ogimet.com/synops.phtml.en

Reforecasts: http://www.esrl.noaa.gov/psd/forecasts/reforecast2/

Hamill TM, Bates GT, Whitaker JS, Murray DR, Fiorino M, Galarneau Jr. TJ, Zhu Y, Lapenta W (2013). NOAA's Second-Generation Global Medium-Range Ensemble Reforecast Data Set. Bulletin of the American Meteorological Society, 94(10), 1553–1565. doi: 10.1175/BAMS-D-12-00014.1

## load data and omit only a couple of observations with some missing values
data("IbkTemperature", package = "lmSubsets")
IbkTemperature <- na.omit(IbkTemperature)

## fit a simple climatological model for the temperature
## with a linear trend and annual/bi-annual harmonic seasonal pattern
m0 <- lm(temp ~ time + sin + cos + sin2 + cos2, data = IbkTemperature)

## fit a simple MOS with 2-meter temperature forecast in addition
## to the climatological model
m1 <- lm(temp ~ t2m + time + sin + cos + sin2 + cos2, data = IbkTemperature)

## graphical comparison and MOS summary
plot(temp ~ time, data = IbkTemperature, type = "l", col = "darkgray")
lines(fitted(m1) ~ time, data = IbkTemperature, col = "darkred")
lines(fitted(m0) ~ time, data = IbkTemperature, lwd = 2)
summary(m1)

## best subset selection of remaining variables for the MOS
## (i.e., forcing the regressors of m1 into the model)
ms2 <- lmSubsets(temp ~ ., data = IbkTemperature,
  include = c("t2m", "time", "sin", "cos", "sin2", "cos2"))
plot(summary(ms2))
image(ms2, size = 8:20)
## -> Note that soil temperature and maximum temperature are selected
## in addition to the 2-meter temperature

## best subset selection of all variables
ms3 <- lmSubsets(temp ~ ., data = IbkTemperature)
plot(summary(ms3))
image(ms3, size = 2:20)
## -> Note that 2-meter temperature is not selected into the best
## BIC model but soil-temperature (and maximum temperature) are used instead

## refit the best BIC subset selections
m2 <- refit(lmSelect(ms2, penalty = "BIC"))
m3 <- refit(lmSelect(ms3, penalty = "BIC"))

## compare BIC
BIC(m0, m1, m2, m3)

## compare RMSE
sqrt(sapply(list(m0, m1, m2, m3), deviance)/nrow(IbkTemperature))

## compare coefficients
cf0 <- coef(m0)
cf1 <- coef(m1)
cf2 <- coef(m2)
cf3 <- coef(m3)
names(cf2) <- gsub("^x", "", names(coef(m2)))
names(cf3) <- gsub("^x", "", names(coef(m3)))
nam <- unique(c(names(cf0), names(cf1), names(cf2), names(cf3)))
cf <- matrix(NA, nrow = length(nam), ncol = 4,
  dimnames = list(nam, c("m0", "m1", "m2", "m3")))
cf[names(cf0), 1] <- cf0
cf[names(cf1), 2] <- cf1
cf[names(cf2), 3] <- cf2
cf[names(cf3), 4] <- cf3
print(round(cf, digits = 3), na.print = "")