.Rhistory

# Perform independent test with known covariance matrix
covMat <- cov(X)  # Example covariance matrix
independent_test(X, covMat = covMa
independent_test <- function(X, covMat = NULL, alpha = NULL) {
independent_test <- function(X, covMat = NULL, alpha = NULL) {
n <- nrow(X)
p <- ncol(X)
# Maximum likelihood estimator of covariance matrix
A <- (n - 1) * cov(X)
sig <- diag(cov(X))
W <- colSums(X)
if (is.null(covMat)) {
# Sigma unknown
test_statistic <- (n - 1) * var(W) / sum(sig)
df <- n - 1
} else {
# Sigma known
test_statistic <- (n - 1) * var(W) / sum(diag(covMat))
df <- n - 1
}
# Define significance levels (alphas)
alphas <- c(0.01, 0.05, 0.01, 0.05, 0.1, 0.2, 0.5)
# Calculate p-values for both tails
p_values_upper <- pchisq(test_statistic, df = df, lower.tail = FALSE)
p_values_lower <- pchisq(test_statistic, df = df, lower.tail = TRUE)
# Calculate final p-value using the minimum of probabilities
p_value <- 2 * pmin(p_values_upper, p_values_lower)
# Create a data frame to store results
results <- data.frame(Test_Statistic = test_statistic,
Alpha = alphas,
P_Value = p_value,
Test_Result = ifelse(p_value <= alphas, "Reject H0", "Fail to Reject H0"))
if (!is.null(alpha)) {
result <- data.frame(Test_Statistic = test_statistic,
Alpha = alpha,
P_Value = p_value,
Test_Result = ifelse(p_value <= alphas, "Reject H0", "Fail to Reject H0"))
return(result)
}
return(results)
}
# Example usage
set.seed(123)  # for reproducibility
X <- matrix(rnorm(100), ncol = 5)  # Example data
# Perform independent test without specifying alpha
independent_test(X)
independent_test <- function(X, covMat = NULL, alpha = NULL) {
n <- nrow(X)
p <- ncol(X)
# Maximum likelihood estimator of covariance matrix
A <- (n - 1) * cov(X)
sig <- diag(cov(X))
W <- colSums(X)
if (is.null(covMat)) {
# Sigma unknown
test_statistic <- (n - 1) * var(W) / sum(sig)
df <- n - 1
} else {
# Sigma known
test_statistic <- (n - 1) * var(W) / sum(diag(covMat))
df <- n - 1
}
# Define significance levels (alphas)
alphas <- c(0.01, 0.05, 0.01, 0.05, 0.1, 0.2, 0.5)
# Calculate p-values for both tails
p_values_upper <- pchisq(test_statistic, df = df, lower.tail = FALSE)
p_values_lower <- pchisq(test_statistic, df = df, lower.tail = TRUE)
# Calculate final p-value using the minimum of probabilities
p_value <- 2 * pmin(p_values_upper, p_values_lower)
# Create a data frame to store results
results <- data.frame(Test_Statistic = test_statistic,
Degrees_of_Freedom = df,
Alpha = alphas,
P_Value = p_value,
Test_Result = ifelse(p_value <= alphas, "Reject H0", "Fail to Reject H0"))
if (!is.null(alpha)) {
result <- data.frame(Test_Statistic = test_statistic,
Degrees_of_Freedom = df,
Alpha = alpha,
P_Value = p_value,
Test_Result = ifelse(p_value <= alphas, "Reject H0", "Fail to Reject H0"))
return(result)
}
return(results)
}
# Example usage
set.seed(123)  # for reproducibility
X <- matrix(rnorm(100), ncol = 5)  # Example data
# Perform independent test without specifying alpha
independent_test(X)
# Perform independent test with alpha = 0.05
independent_test(X, alpha = 0.05)
covMat <- covMatCS(5)
covMat
covMatC <- function(p, sigma2 = 1, rho = NULL) {
if (is.null(rho)) {
rho <- stats::runif(1) # If rho is not provided, generate a random value between 0 and 1
}
# Create an empty matrix
matrix <- matrix(0, nrow = p, ncol = p)
# Populate the covariance matrix
for (i in 1:p) {
for (j in 1:p) {
if (i == j) {
matrix[i, j] <- sigma2
} else {
distance <- min(abs(i - j), p - abs(i - j))
matrix[i, j] <- sigma2 * exp(-rho * distance)
}
}
}
if (matrixcalc::is.positive.definite(matrix)) {
message <- "The covariance matrix is positive definite."
} else {
message <- "The covariance matrix is not positive definite."
}
cat(message, "\n")
return(matrix)
}
covMatC(p = 5)
covMatC <- function(p, sigma2 = 1, rho = NULL) {
if (is.null(rho)) {
rho <- stats::runif(1) # If rho is not provided, generate a random value between 0 and 1
}
# Create an empty matrix
matrix <- matrix(0, nrow = p, ncol = p)
# Populate the covariance matrix
for (i in 1:p) {
for (j in 1:p) {
if (i == j) {
matrix[i, j] <- sigma2
} else {
distance <- min(abs(i - j), p - abs(i - j))
matrix[i, j] <- sigma2 * exp(-rho * distance)
}
}
}
if (matrixcalc::is.positive.definite(matrix)) {
message <- "The covariance matrix is positive definite."
} else {
message <- "The covariance matrix is not positive definite."
}
cat(message, "\n")
cat("\n")
return(matrix)
}
covMatC(5)
mat = covMatC(5)
mat
covMat <- covMatCS(5)
covMat
covMatCS(5)
covMat <- covMatC(5)
covMat
covMat <- covMatC(5)
covMat
X <- mvrnorm(n = 100, mu = rep(0,5), Sigma = covMat)
X <- mvrnorm(n = 100, mu = rep(0,5), Sigma = covMat)
X <- MASS::mvrnorm(n = 100, mu = rep(0,5), Sigma = covMat)
# Perform independent test with alpha = 0.05
independent_test(X, alpha = 0.05)
covMat <- diag(5)
X <- MASS::mvrnorm(n = 100, mu = rep(0,5), Sigma = covMat)
# Perform independent test with alpha = 0.05
independent_test(X, alpha = 0.05)
library(indepTest)
library(indepTest)
library(indepTest)
library(indepTest)
library(indepTest)
cov_mat
library(MASS)
n = 100 # Sample Size
p = 5
rho = 0.014
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
cov_mat
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
#'
#' #' # Covariance structure with Circular structure
#' cov_mat <- covMatC(p = p, rho = rho)
#' data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
#' indTest(data)
#' indTest(data, covMat = cov_mat)
#'
#'
#' @export
#' @importFrom stats cov var pchisq
indTest <- function(X, covMat = NULL, alpha = 0.05) {
n <- nrow(X)
p <- ncol(X)
# Maximum likelihood estimator of covariance matrix
A <- (n - 1) * stats::cov(X)
if (is.null(covMat)) {
covMat = stats::cov(X)
}
off_diag <- covMat[!row(covMat) == col(covMat)]
if (any(off_diag > 0) && any(off_diag < 0)) {
cat("Alert:")
cat('\n')
cat("The off-diagonal elements of the covariance matrix are not the same signs. The results may not be reliable.")
cat("\n")
}
# Check if the covariance matrix is symmetric
if (!isSymmetric(covMat)) {
stop("Covariance matrix is not symmetric.")
}
sig <- diag(cov(X))
W <- rowSums(X)
if (is.null(covMat)) {
# Sigma unknown
test_statistic <- (n - 1) * stats::var(W) / sum(sig)
df <- n - 1
} else {
# Sigma known
test_statistic <- (n - 1) * stats::var(W) / sum(diag(covMat))
df <- n - 1
}
# Calculate p-values for both tails
p_values_upper <- stats::pchisq(test_statistic, df = df, lower.tail = FALSE)
p_values_lower <- stats::pchisq(test_statistic, df = df, lower.tail = TRUE)
# Calculate final p-value using the minimum of probabilities
p_value <- 2 * min(p_values_upper, p_values_lower)
# Create a data frame to store results
result <- data.frame(Test_Statistic = test_statistic,
Degrees_of_Freedom = df,
Alpha = alpha,
P_Value = p_value,
Test_Result = ifelse(p_value <= alpha, "Reject H0", "Fail to Reject H0"))
return(result)
}
indTest(data)
indTest(data, covMat = cov_mat)
#' # Covariance structure with Circular structure
cov_mat <- covMatC(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
indTest(data, covMat = cov_mat)
cov_mat
n = 10 # Sample Size
p = 20  # number of variables
rho = 0.014
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
data
indTest(data)
rho = 0.4
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
indTest(data)
cov_mat)
cov_mat
rho = 0.9
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
cov_mat
rho = 0.4
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
cov_mat
cov_mat
cov_mat <- covMatAR(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
indTest(data)
cov(data)
indTest(data, covMat = cov_mat)
usethis::use_git_config(user.name = "YourName", user.email = "your@mail.com")
install.packages("usethis")
library(usethis)
usethis::use_git_config(user.name = "fjmmarques", user.email = "fjm@fct.unl.pt")
n = 50 # Sample Size
p = 5  # number of variables
rho = 0.4
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
indTest(data, covMat = cov_mat)
# data with missing vales
missing_rate <- 0.1  # 10% missing values
missing_index_row <- sample(1:n, size = round(n * missing_rate))
missing_index_col <- sample(1:p, size = 1)
data[missing_index_row, missing_index_col] <- NA  # Introduce missing values
indTest(data)
indTest(data, covMat = cov_mat)
# data with missing vales
missing_rate <- 0.1  # 10% missing values
missing_index_row <- sample(1:n, size = round(n * missing_rate))
missing_index_col <- sample(1:p, size = 1)
data[missing_index_row, missing_index_col] <- NA  # Introduce missing values
indTest(data)
#'
#' #' # Covariance structure with Circular structure
#' cov_mat <- covMatC(p = p, rho = rho)
#' data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
#' indTest(data)
#' indTest(data, covMat = cov_mat)
#'
#'
#' @export
#' @importFrom stats cov var pchisq
indTest <- function(X, covMat = NULL, alpha = 0.05) {
if (any(is.na(X))) {
cat("Alert:")
cat('\n')
cat("The data has missing values. The missing values are handled by casewise deletion (and if there are no complete cases, that gives an error)")
cat("\n")
# [We should discuss the following]
# cov(X, use = 'complete.ons')
X <- data.frame(na.omit(X))
}
n <- nrow(X)
p <- ncol(X)
# Maximum likelihood estimator of covariance matrix
A <- (n - 1) * stats::cov(X)
if (is.null(covMat)) {
covMat = stats::cov(X)
}
off_diag <- covMat[!row(covMat) == col(covMat)]
if (any(off_diag > 0) && any(off_diag < 0)) {
cat("Alert:")
cat('\n')
cat("The off-diagonal elements of the covariance matrix are not the same signs. The results may not be reliable.")
cat("\n")
}
# Check if the covariance matrix is symmetric
if (!isSymmetric(covMat)) {
stop("Covariance matrix is not symmetric.")
}
sig <- diag(cov(X))
W <- rowSums(X)
if (is.null(covMat)) {
# Sigma unknown
test_statistic <- (n - 1) * stats::var(W) / sum(sig)
df <- n - 1
} else {
# Sigma known
test_statistic <- (n - 1) * stats::var(W) / sum(diag(covMat))
df <- n - 1
}
# Calculate p-values for both tails
p_values_upper <- stats::pchisq(test_statistic, df = df, lower.tail = FALSE)
p_values_lower <- stats::pchisq(test_statistic, df = df, lower.tail = TRUE)
# Calculate final p-value using the minimum of probabilities
p_value <- 2 * min(p_values_upper, p_values_lower)
# Create a data frame to store results
result <- data.frame(Test_Statistic = test_statistic,
Degrees_of_Freedom = df,
Alpha = alpha,
P_Value = p_value,
Test_Result = ifelse(p_value <= alpha, "Reject H0", "Fail to Reject H0"))
return(result)
}
#'
#' #' # Covariance structure with Circular structure
#' cov_mat <- covMatC(p = p, rho = rho)
#' data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
#' indTest(data)
#' indTest(data, covMat = cov_mat)
#'
#'
#' @export
#' @importFrom stats cov var pchisq
indTest <- function(X, covMat = NULL, alpha = 0.05) {
if (any(is.na(X))) {
cat("Alert:")
cat('\n')
cat("The data has missing values. The missing values are handled by casewise deletion (and if there are no complete cases, that gives an error)")
cat("\n")
# [We should discuss the following]
# cov(X, use = 'complete.ons')
X <- data.frame(na.omit(X))
}
n <- nrow(X)
p <- ncol(X)
# Maximum likelihood estimator of covariance matrix
A <- (n - 1) * stats::cov(X)
if (is.null(covMat)) {
covMat = stats::cov(X)
}
off_diag <- covMat[!row(covMat) == col(covMat)]
if (any(off_diag > 0) && any(off_diag < 0)) {
cat("Alert:")
cat('\n')
cat("The off-diagonal elements of the covariance matrix are not the same signs. The results may not be reliable.")
cat("\n")
}
# Check if the covariance matrix is symmetric
if (!isSymmetric(covMat)) {
stop("Covariance matrix is not symmetric.")
}
sig <- diag(cov(X))
W <- rowSums(X)
if (is.null(covMat)) {
# Sigma unknown
test_statistic <- (n - 1) * stats::var(W) / sum(sig)
df <- n - 1
} else {
# Sigma known
test_statistic <- (n - 1) * stats::var(W) / sum(diag(covMat))
df <- n - 1
}
# Calculate p-values for both tails
p_values_upper <- stats::pchisq(test_statistic, df = df, lower.tail = FALSE)
p_values_lower <- stats::pchisq(test_statistic, df = df, lower.tail = TRUE)
# Calculate final p-value using the minimum of probabilities
p_value <- 2 * min(p_values_upper, p_values_lower)
# Create a data frame to store results
result <- data.frame(Test_Statistic = test_statistic,
Degrees_of_Freedom = df,
Alpha = alpha,
P_Value = p_value,
Test_Result = ifelse(p_value <= alpha, "Reject H0", "Fail to Reject H0"))
return(result)
}
library(MASS)
n = 50 # Sample Size
p = 5  # number of variables
rho = 0.4
# Covariance structure with Autoregressive structure
cov_mat <- covMatAR(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
indTest(data, covMat = cov_mat)
# data with missing vales
missing_rate <- 0.1  # 10% missing values
missing_index_row <- sample(1:n, size = round(n * missing_rate))
missing_index_col <- sample(1:p, size = 1)
data[missing_index_row, missing_index_col] <- NA  # Introduce missing values
indTest(data)
indTest(data, covMat = cov_mat)
# Covariance structure with Compound Symmetry structure
cov_mat <- covMatCS(p = p, rho = rho)
data <- mvrnorm(n = n, mu = rep(0,p), Sigma = cov_mat)
indTest(data)
indTest(data, covMat = cov_mat)
# install.packages("ggplot2")
# install.packages("maps")
library(ggplot2)
library(maps)
# Import the data with coordinates
world <- map_data("world")
# Plot the map. group = group connects the points in the correct order
ggplot(data = world, aes(x = long, y = lat, group = group)) +
geom_polygon()
# Import the data with coordinates
world <- map_data("portugal")
# Plot the map. group = group connects the points in the correct order
ggplot(data = port, aes(x = long, y = lat, group = group)) +
geom_polygon()
# Import the data with coordinates
port <- map_data("portugal")
# Import the data with coordinates
port <- map_data("port")
# install.packages("ggplot2")
# install.packages("maps")
library(ggplot2)
library(maps)
#Add the data to the shapefile
europe_df <- europe %>% left_join(df, by=c("name_long" ="Country"))
#Add the data to the shapefile
europe_df <- europe %>% left_join(df, by=c("name_long" ="Portugal"))
#Add the data to the shapefile
europe_df <-left_join(df, by=c("name_long" ="Portugal"))
library(maps)
library(ggplot2)
## make a df with only the country to overlap
map_data_es <- map_data('world')[map_data('world')$region == "Spain",]
## The map (maps + ggplot2 )
ggplot() +
## First layer: worldwide map
geom_polygon(data = map_data("world"),
aes(x=long, y=lat, group = group),
color = '#9c9c9c', fill = '#f3f3f3') +
## Second layer: Country map
geom_polygon(data = map_data_es,
aes(x=long, y=lat, group = group),
color = 'red', fill = 'pink') +
coord_map() +
coord_fixed(1.3,
xlim = c(-13.5, 8.5),
ylim = c(34, 45)) +
ggtitle("A map of Spain") +
theme(panel.background =element_rect(fill = 'blue'))
## make a df with only the country to overlap
map_data_es <- map_data('world')[map_data('world')$region == "Portugal",]
## make a df with only the country to overlap
map_data_es <- map_data('world')[map_data('world')$region == "Portugal",]
## The map (maps + ggplot2 )
ggplot() +
## First layer: worldwide map
geom_polygon(data = map_data("world"),
aes(x=long, y=lat, group = group),
color = '#9c9c9c', fill = '#f3f3f3') +
## Second layer: Country map
geom_polygon(data = map_data_es,
aes(x=long, y=lat, group = group),
color = 'red', fill = 'pink') +
coord_map() +
coord_fixed(1.3,
xlim = c(-13.5, 8.5),
ylim = c(34, 45)) +
ggtitle("A map of Spain") +
theme(panel.background =element_rect(fill = 'blue'))
## make a df with only the country to overlap
map_data_es <- map_data('portugal')[map_data('world')$region == "Portugal",]