# Some of the examples may run some time.
# Canadian weather data set
# There are three samples of mean temperature and precipitations for
# fifteen weather stations in Western Canada, another fifteen in Eastern Canada,
# and the remaining five in Northern Canada.
# one functional variable - temperature
library(fda)
data_set_t <- t(CanadianWeather$dailyAv[,, "Temperature.C"])
# number of samples
k <- 3
# number of variables
p <- 1
# preparing data set
gr_label <- rep(c(1, 2, 3), c(15, 15, 5))
data_set <- list(list(data_set_t[gr_label == 1, ]),
list(data_set_t[gr_label == 2, ]),
list(data_set_t[gr_label == 3, ]))
# Tukey's contrast matrix
h_tukey <- GFDmcv::contr_mat(k, type = "Tukey")
h_tukey_m <- kronecker(h_tukey, diag(p))
# plots for pointwise Hotelling's T^2-test statistics
oldpar <- par(mfrow = c(2, 2), mar = c(4, 2, 2, 0.1))
plot(ph_test_statistic(data_set, h_tukey_m), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Global hypothesis", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_tukey_m[1, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Contrast 1", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_tukey_m[2, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Contrast 2", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_tukey_m[3, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Contrast 3", xlab = "Day")
par(oldpar)
# Dunnett's contrast matrix
h_dunnett <- GFDmcv::contr_mat(k, type = "Dunnett")
h_dunnett_m <- kronecker(h_dunnett, diag(p))
# plots for pointwise Hotelling's T^2-test statistics
oldpar <- par(mfrow = c(3, 1), mar = c(4, 2, 2, 0.1))
plot(ph_test_statistic(data_set, h_dunnett_m), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
main = "Global hypothesis", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_dunnett_m[1, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
main = "Contrast 1", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_dunnett_m[2, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
main = "Contrast 2", xlab = "Day")
par(oldpar)
# two functional variables - temperature and precipitation
library(fda)
data_set_t <- t(CanadianWeather$dailyAv[,, "Temperature.C"])
data_set_p <- t(CanadianWeather$dailyAv[,, "Precipitation.mm"])
# number of samples
k <- 3
# number of variables
p <- 2
# preparing data set
gr_label <- rep(c(1, 2, 3), c(15, 15, 5))
data_set <- list(list(data_set_t[gr_label == 1, ], data_set_p[gr_label == 1, ]),
list(data_set_t[gr_label == 2, ], data_set_p[gr_label == 2, ]),
list(data_set_t[gr_label == 3, ], data_set_p[gr_label == 3, ]))
# Tukey's contrast matrix
h_tukey <- GFDmcv::contr_mat(k, type = "Tukey")
h_tukey_m <- kronecker(h_tukey, diag(p))
# plots for pointwise Hotelling's T^2-test statistics
oldpar <- par(mfrow = c(2, 2), mar = c(4, 2, 2, 0.1))
plot(ph_test_statistic(data_set, h_tukey_m), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Global hypothesis", xlab = "Day")
plot(ph_test_statistic(data_set, h_tukey_m[1:2, ]), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Contrast 1", xlab = "Day")
plot(ph_test_statistic(data_set, h_tukey_m[3:4, ]), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Contrast 2", xlab = "Day")
plot(ph_test_statistic(data_set, h_tukey_m[5:6, ]), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
main = "Contrast 3", xlab = "Day")
par(oldpar)
# Dunnett's contrast matrix
h_dunnett <- GFDmcv::contr_mat(k, type = "Dunnett")
h_dunnett_m <- kronecker(h_dunnett, diag(p))
# plots for pointwise Hotelling's T^2-test statistics
oldpar <- par(mfrow = c(3, 1), mar = c(4, 2, 2, 0.1))
plot(ph_test_statistic(data_set, h_dunnett_m), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
main = "Global hypothesis", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_dunnett_m[1, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
main = "Contrast 1", xlab = "Day")
plot(ph_test_statistic(data_set, matrix(h_dunnett_m[2, ], 1)), type = "l",
ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
main = "Contrast 2", xlab = "Day")
par(oldpar)
Run the code above in your browser using DataLab