beav2.Rd
Reynolds (1994) describes a small part of a study of the long-term temperature dynamics of beaver Castor canadensis in north-central Wisconsin. Body temperature was measured by telemetry every 10 minutes for four females, but data from a one period of less than a day for each of two animals is used there.
beav2
The beav2
data frame has 100 rows and 4 columns.
This data frame contains the following columns:
day
Day of observation (in days since the beginning of 1990), November 3–4.
time
Time of observation, in the form 0330
for 3.30am.
temp
Measured body temperature in degrees Celsius.
activ
Indicator of activity outside the retreat.
P. S. Reynolds (1994) Time-series analyses of beaver body temperatures. Chapter 11 of Lange, N., Ryan, L., Billard, L., Brillinger, D., Conquest, L. and Greenhouse, J. eds (1994) Case Studies in Biometry. New York: John Wiley and Sons.
Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer.
attach(beav2)
beav2$hours <- 24*(day-307) + trunc(time/100) + (time%%100)/60
plot(beav2$hours, beav2$temp, type = "l", xlab = "time",
ylab = "temperature", main = "Beaver 2")
usr <- par("usr"); usr[3:4] <- c(-0.2, 8); par(usr = usr)
lines(beav2$hours, beav2$activ, type = "s", lty = 2)
temp <- ts(temp, start = 8+2/3, frequency = 6)
activ <- ts(activ, start = 8+2/3, frequency = 6)
acf(temp[activ == 0]); acf(temp[activ == 1]) # also look at PACFs
ar(temp[activ == 0]); ar(temp[activ == 1])
#>
#> Call:
#> ar(x = temp[activ == 0])
#>
#> Coefficients:
#> 1
#> 0.7392
#>
#> Order selected 1 sigma^2 estimated as 0.02011
#>
#> Call:
#> ar(x = temp[activ == 1])
#>
#> Coefficients:
#> 1
#> 0.7894
#>
#> Order selected 1 sigma^2 estimated as 0.01792
arima(temp, order = c(1,0,0), xreg = activ)
#>
#> Call:
#> arima(x = temp, order = c(1, 0, 0), xreg = activ)
#>
#> Coefficients:
#> ar1 intercept activ
#> 0.8733 37.1920 0.6139
#> s.e. 0.0684 0.1187 0.1381
#>
#> sigma^2 estimated as 0.01518: log likelihood = 66.78, aic = -125.55
dreg <- cbind(sin = sin(2*pi*beav2$hours/24), cos = cos(2*pi*beav2$hours/24))
arima(temp, order = c(1,0,0), xreg = cbind(active=activ, dreg))
#>
#> Call:
#> arima(x = temp, order = c(1, 0, 0), xreg = cbind(active = activ, dreg))
#>
#> Coefficients:
#> ar1 intercept active dreg.sin dreg.cos
#> 0.7905 37.1674 0.5322 -0.282 0.1201
#> s.e. 0.0681 0.0939 0.1282 0.105 0.0997
#>
#> sigma^2 estimated as 0.01434: log likelihood = 69.83, aic = -127.67
## IGNORE_RDIFF_BEGIN
library(nlme) # for gls and corAR1
beav2.gls <- gls(temp ~ activ, data = beav2, correlation = corAR1(0.8),
method = "ML")
summary(beav2.gls)
#> Generalized least squares fit by maximum likelihood
#> Model: temp ~ activ
#> Data: beav2
#> AIC BIC logLik
#> -125.5505 -115.1298 66.77523
#>
#> Correlation Structure: AR(1)
#> Formula: ~1
#> Parameter estimate(s):
#> Phi
#> 0.8731771
#>
#> Coefficients:
#> Value Std.Error t-value p-value
#> (Intercept) 37.19195 0.1131328 328.7460 0
#> activ 0.61418 0.1087286 5.6487 0
#>
#> Correlation:
#> (Intr)
#> activ -0.582
#>
#> Standardized residuals:
#> Min Q1 Med Q3 Max
#> -2.42080776 -0.61510520 -0.03573836 0.81641138 2.15153496
#>
#> Residual standard error: 0.2527856
#> Degrees of freedom: 100 total; 98 residual
summary(update(beav2.gls, subset = 6:100))
#> Generalized least squares fit by maximum likelihood
#> Model: temp ~ activ
#> Data: beav2
#> Subset: 6:100
#> AIC BIC logLik
#> -124.981 -114.7654 66.49048
#>
#> Correlation Structure: AR(1)
#> Formula: ~1
#> Parameter estimate(s):
#> Phi
#> 0.8380448
#>
#> Coefficients:
#> Value Std.Error t-value p-value
#> (Intercept) 37.25001 0.09634047 386.6496 0
#> activ 0.60277 0.09931904 6.0690 0
#>
#> Correlation:
#> (Intr)
#> activ -0.657
#>
#> Standardized residuals:
#> Min Q1 Med Q3 Max
#> -2.0231494 -0.8910348 -0.1497564 0.7640939 2.2719468
#>
#> Residual standard error: 0.2188542
#> Degrees of freedom: 95 total; 93 residual
detach("beav2"); rm(temp, activ)
## IGNORE_RDIFF_END