Model robust and potential model misspecification for Electric Consumption Data
Source:vignettes/MR_MM_Linear_Reg_Electricity_Con.Rmd
MR_MM_Linear_Reg_Electricity_Con.Rmd
# Load the R packages
library(NeEDS4BigData)## Error in get(paste0(generic, ".", class), envir = get_method_env()) :
## object 'type_sum.accel' not found
library(ggplot2)
library(ggpubr)
library(kableExtra)
library(tidyr)
library(psych)
library(gam)
# Theme for plots
Theme_special<-function(){
theme(legend.key.width= unit(1,"cm"),
axis.text.x= element_text(color= "black", size= 12, angle= 30, hjust=0.75),
axis.text.y= element_text(color= "black", size= 12),
strip.text= element_text(colour= "black", size= 12, face= "bold"),
panel.grid.minor.x= element_blank(),
axis.title= element_text(color= "black", face= "bold", size= 12),
legend.text= element_text(color= "black", size= 11),
legend.position= "bottom",
legend.margin= margin(0,0,0,0),
legend.box.margin= margin(-1,-2,-1,-2))
}Understanding the electric consumption data
``Electric power consumption’’ data (Hebrail and Berard 2012), which contains \(2,049,280\) measurements for a house located at Sceaux, France between December 2006 and November 2010. The data contains \(4\) columns and has \(2,049,280\) observations, first column is the response variable and the rest are covariates, however we only use the first \(10\%\) of the data for this explanation. The response \(y\) is the log scaled intensity, while the covariates are active electrical energy in the a) voltage (\(X_1\)) b) kitchen (\(X_2\)), c) laundry room (\(X_3\)) and d) water-heater and air-conditioner (\(X_4\)). The covariates are scaled to be mean of zero and variance of one.
Given data is analysed under two different scenarios,
- model robust or average subsampling methods assuming that a set of models can describe the data.
- subsampling method assuming the main effects model is potentially misspecified.
# Selecting 25% of the big data and prepare it
indexes<-1:ceiling(nrow(Electric_consumption)*0.10)
Original_Data<-cbind(Electric_consumption[indexes,1],1,
Electric_consumption[indexes,-1])
colnames(Original_Data)<-c("Y",paste0("X",0:ncol(Original_Data[,-c(1,2)])))
# Scaling the covariate data
for (j in 3:6) {
Original_Data[,j]<-scale(Original_Data[,j])
}
head(Electric_consumption) %>%
kable(format = "html",
caption = "First few observations of the electric consumption data.")| Intensity | Voltage | EE_Kitchen | EE_Laundry | EE_WH_AC |
|---|---|---|---|---|
| 2.912351 | 234.84 | 0 | 1 | 17 |
| 3.135494 | 233.63 | 0 | 1 | 16 |
| 3.135494 | 233.29 | 0 | 2 | 17 |
| 3.135494 | 233.74 | 0 | 1 | 17 |
| 2.760010 | 235.68 | 0 | 1 | 17 |
| 2.708050 | 235.02 | 0 | 2 | 17 |
# Setting the sample sizes
N<-nrow(Original_Data); M<-250; k<-seq(8,20,by=2)*100; rep_k<-rep(k,each=M)Model robust or average subsampling
The method \(A\)- and \(L\)-optimality of model robust or average
subsampling (Mahendran, Thompson, and McGree
2023) is compared against the \(A\)- and \(L\)-optimality subsampling (Ai et al. 2021; Yao and Wang 2021) method. Here
five models are selected out of sixteen models, which are 1) main
effects model (\(\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+\beta_4X_4\)),
2-5) main effects model with the squared term of each covariate (\(X^2_1 / X^2_2 / X^2_3 / X^2_4\)), 6-11)
main effects model with two squared terms of the covariates, 12-15) main
effects model with three squared terms of the covariates, and 16) main
effects model with all squared terms. The five models are selected based
on the smallest AIC values.
For each model \(l\) the mean squared
error of the model parameters \(MSE_l(\tilde{\beta}_k,\hat{\beta})=\frac{1}{MJ}
\sum_{i=1}^M \sum_{j=1}^J (\tilde{\beta}_{k,j} -
\hat{\beta}_j)^2\) are calculated for the \(M=250\) simulations across the sample sizes
\(k=r_f=(800,\ldots,2000)\) and the
initial sample size is \(r_0=400\).
Here, for the \(l\)-th model \(\tilde{\beta}_k\) is the estimated model
parameters from the sample of size \(k\) and \(\hat{\beta}\) is the estimated model
parameters from the full big data, while \(j\) is index of the model parameter.
# Define the subsampling methods and their respective colours, shapes and line types
Method_Names<-c("A-Optimality","L-Optimality","A-Optimality MR","L-Optimality MR")
Method_Colour<-c("#D82632","#A50021","#50FF50","#00BB00")
Method_Shape_Types<-c(rep(8,2),rep(17,2))
# Preparing the data for the model average method with squared terms
No_of_Variables<-ncol(Original_Data[,-c(1,2)])
Squared_Terms<-paste0("X",1:No_of_Variables,"^2")
term_no <- 2
All_Models <- list(c("X0",paste0("X",1:No_of_Variables)))
Original_Data_ModelRobust<-cbind(Original_Data,Original_Data[,-c(1,2)]^2)
colnames(Original_Data_ModelRobust)<-c("Y","X0",paste0("X",1:No_of_Variables),
paste0("X",1:No_of_Variables,"^2"))
Temp_Data<-Original_Data_ModelRobust[sample(1:N,1000),]
AIC_Values<-NULL
model <- lm(Y~.-1, data = Temp_Data[,c("Y","X0",paste0("X",1:No_of_Variables))])
AIC_Values[1]<-AIC(model)
for (i in 1:No_of_Variables)
{
x <- as.vector(combn(Squared_Terms,i,simplify = FALSE))
for(j in 1:length(x))
{
All_Models[[term_no]] <- c("X0",paste0("X",1:No_of_Variables),x[[j]])
temp_data<-as.data.frame(Temp_Data[,c("Y",All_Models[[term_no]])])
model <- lm(Y~.-1, data = temp_data)
AIC_Values[term_no]<-AIC(model)
term_no <- term_no+1
}
}
Model_Set<-5
Best_Indices <- order(AIC_Values,decreasing = FALSE)[1:Model_Set]
Model_Apriori<-rep(1/Model_Set,Model_Set)
All_Models<-All_Models[Best_Indices]
names(All_Models)<-paste0("Model_",1:length(All_Models))Covariates in the selected five models
- X0, X1, X2, X3, X4, X1^2, X3^2, X4^2
- X0, X1, X2, X3, X4, X3^2, X4^2
- X0, X1, X2, X3, X4, X3^2
- X0, X1, X2, X3, X4, X1^2, X3^2
- X0, X1, X2, X3, X4, X1^2, X2^2, X3^2, X4^2
A priori probabilities are equal
Consider for \(Q=5\) each model has an equal a priori probability (i.e \(\alpha_q=1/5,q=1,\ldots,5\)). Below is the code of implementation for this scenario.
All_Covariates<-colnames(Original_Data_ModelRobust)[-1]
# A- and L-optimality model robust subsampling for linear regression
NeEDS4BigData::modelRobustLinSub(r0=400,rf=rep_k,
Y=as.matrix(Original_Data_ModelRobust[,1]),
X=as.matrix(Original_Data_ModelRobust[,-1]),
N=N,Apriori_probs=rep(1/Model_Set,Model_Set),
All_Combinations=All_Models,
All_Covariates=All_Covariates)->Results## Step 1 of the algorithm completed.## Step 2 of the algorithm completed.
Final_Beta_modelRobust<-Results$Beta_Estimates
# Mean squared error and their respective plots for all five models
plot_list_MR<-list()
for (i in 1:length(All_Models)) {
lm(Y~.-1,data=Original_Data_ModelRobust[,c("Y",All_Models[[i]])])->All_Results
All_Beta<-coefficients(All_Results)
matrix(rep(All_Beta,by=nrow(Final_Beta_modelRobust[[i]])),
nrow = nrow(Final_Beta_modelRobust[[i]]),
ncol = ncol(Final_Beta_modelRobust[[i]][,-c(1,2)]),byrow = TRUE)->All_Beta
data.frame("Method"=Final_Beta_modelRobust[[i]]$Method,
"Sample"=Final_Beta_modelRobust[[i]]$rf,
"SE"=rowSums((All_Beta-
Final_Beta_modelRobust[[i]][,-c(1,2)])^2))->MSE_Beta_MR
Mean_Data_MR <- MSE_Beta_MR |>
group_by(Method,Sample) |>
dplyr::summarise(Mean = mean(SE),
min=quantile(SE,0.05),
max=quantile(SE,0.95), .groups = "drop")
ggplot(Mean_Data_MR,aes(x=factor(Sample),y=Mean,color=Method,shape=Method))+
xlab("Sample size")+ylab("Squared Error")+
geom_point(size=3,position = position_dodge(width = 0.5))+
geom_line(data=Mean_Data_MR,aes(x=factor(Sample),y=Mean,group=Method,color=Method),
position = position_dodge(width = 0.5))+
geom_errorbar(data=Mean_Data_MR,aes(ymin=min,ymax=max),
width=0.3,position = position_dodge(width = 0.5))+
ggtitle(paste0("Model ",i))+
scale_color_manual(values = Method_Colour)+
scale_shape_manual(values= Method_Shape_Types)+
theme_bw()+guides(colour= guide_legend(nrow = 2))+
Theme_special()->plot_list_MR[[i]]
}
ggarrange(plotlist = plot_list_MR,nrow = 3,ncol = 2,labels = "auto",
common.legend = TRUE,legend = "bottom")
Mean squared error for all the models with equal apriori in the order a to e for Model 1 to 5 across the subsampling methods under comparison.
Main effects model is potentially misspecified
The final scenario is for comparison of the subsampling method under the assumption that the main effects model is potentially misspecified against the \(A\)- and \(L\)-optimality subsampling method. Under the subsampling method that accounts for potential model misspecification we take the scaling factor of \(\alpha=10\). As in the above scenarios the number of simulations and the sample sizes stay the same. We compare the mean squared error of the estimated model parameters, however as we assume the model is potentially misspecified the asymptotic approximation of the mean squared error from the predictions are calculated as well. Below is the code for this implementation.
# Define the subsampling methods and their respective colours, shapes and line types
Method_Names<-c("A-Optimality","L-Optimality","L1-Optimality","RLmAMSE",
"RLmAMSE Log Odds 10","RLmAMSE Power 10")
Method_Colour<-c("#D82632","#A50021","greenyellow","#BBFFBB","#50FF50","#00BB00")
Method_Shape_Types<-c(rep(8,3),rep(17,3))
# A- and L-optimality and RLmAMSE model misspecified subsampling for linear regression
NeEDS4BigData::modelMissLinSub(r0=400,rf=rep_k,
Y=as.matrix(Original_Data[,1]),
X=as.matrix(Original_Data[,-1]),
N=N,Alpha=10,proportion=0.5)->Results## 50% or >=50% of the big data is used to help find AMSE for the subsamples,
## this could take some time.## Warning: executing %dopar% sequentially: no parallel backend registered## Step 1 of the algorithm completed.## Step 2 of the algorithm completed.
Final_Beta_modelMiss<-Results$Beta_Estimates
Final_AMSE_modelMiss<-Results$AMSE_Estimates
lm(Y~.-1,data = Original_Data)->All_Results
coefficients(All_Results)->All_Beta
matrix(rep(All_Beta,by=nrow(Final_Beta_modelMiss)),nrow = nrow(Final_Beta_modelMiss),
ncol = ncol(Final_Beta_modelMiss[,-c(1,2)]),byrow = TRUE)->All_Beta
# Plots for the mean squared error of the model parameter estimates
# and the AMSE for the main effects model
MSE_Beta_modelMiss<-data.frame("Method"=Final_Beta_modelMiss$Method,
"Sample"=Final_Beta_modelMiss$rf,
"SE"=rowSums((All_Beta -
Final_Beta_modelMiss[,-c(1,2)])^2))
Mean_Data_MM <- MSE_Beta_modelMiss |>
group_by(Method,Sample) |>
dplyr::summarise(Mean = mean(SE),
min=quantile(SE,0.05),
max=quantile(SE,0.95), .groups = "drop")
ggplot(Mean_Data_MM,aes(x=factor(Sample),y=Mean,color=Method,shape=Method))+
xlab("Sample size")+ylab("MSE")+
geom_point(size=3,position = position_dodge(width = 0.5))+
geom_line(data=Mean_Data_MM,aes(x=factor(Sample),y=Mean,group=Method,color=Method),
position = position_dodge(width = 0.5))+
geom_errorbar(data=Mean_Data_MM,aes(ymin=min,ymax=max),
width=0.3,position = position_dodge(width = 0.5))+
scale_color_manual(values = Method_Colour)+
scale_shape_manual(values = Method_Shape_Types)+
theme_bw()+guides(colour = guide_legend(nrow = 3))+Theme_special()->p1
Mean_AMSE<-Final_AMSE_modelMiss |> dplyr::group_by(rf,Method) |>
dplyr::summarise(Mean=mean(AMSE),
min=quantile(AMSE,0.05),
max=quantile(AMSE,0.95),.groups ='drop')
ggplot(Mean_AMSE,aes(x=factor(rf),y=Mean,color=Method,shape=Method)) +
geom_point(size=3,position = position_dodge(width = 0.5))+
geom_line(data=Mean_AMSE,aes(x=factor(rf),y=Mean,group=Method,color=Method),
position = position_dodge(width = 0.5))+
geom_errorbar(data=Mean_AMSE,aes(ymin=min,ymax=max),
width=0.3,position = position_dodge(width = 0.5))+
xlab("Sample size")+ylab("Mean AMSE")+
scale_color_manual(values = Method_Colour)+
scale_shape_manual(values = Method_Shape_Types)+
theme_bw()+guides(colour = guide_legend(nrow = 3))+Theme_special()->p2
ggarrange(p1,p2,nrow = 2,ncol = 1,labels = "auto")
MSE for model parameters (a) and AMSE (b) for the potentially misspecified main effects model across the subsampling methods under comparison, with 5% and 95% percentile intervals.