First example
Aim
We want to creates \(NDatasets=200\) datasets with \(\textrm{length}(group)=10\) variables and \(N=10\) observations. In that example we want \(9\) groups:
- \(x_1\) and \(x_{10}\) belong to the first group and the intra-group Pearson correlation for this group is equal to \(.95\),
- \(x_2\) belongs to the second group,
- \(x_3\) belongs to the third group,
- …
- \(x_9\) belongs to the ninth group.
Correlation structure
The correlation structure of the explanatory variables of the dataset
is provided by group and the intra-group Pearson
correlation value for each of the groups by cor_group. A
value must be provided even for single variable groups and the number of
variables is length of the group vector. Use the
simulation_cor function to create the correlation matrix
(CM).
library(SelectBoost)
group<-c(1:9,1) #10 variables
cor_group<-rep(0.95,9)
CM<-simulation_cor(group,cor_group)
CM
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 1.00 0 0 0 0 0 0 0 0 0.95
#> [2,] 0.00 1 0 0 0 0 0 0 0 0.00
#> [3,] 0.00 0 1 0 0 0 0 0 0 0.00
#> [4,] 0.00 0 0 1 0 0 0 0 0 0.00
#> [5,] 0.00 0 0 0 1 0 0 0 0 0.00
#> [6,] 0.00 0 0 0 0 1 0 0 0 0.00
#> [7,] 0.00 0 0 0 0 0 1 0 0 0.00
#> [8,] 0.00 0 0 0 0 0 0 1 0 0.00
#> [9,] 0.00 0 0 0 0 0 0 0 1 0.00
#> [10,] 0.95 0 0 0 0 0 0 0 0 1.00Explanatory dataset generation
Then generation of an explanatory dataset with \(N=10\) observations is made by the
simulation_X function.
set.seed(3141)
N<-10
X<-simulation_X(N,CM)
X
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.69029599 -1.7343187 0.38993973 0.7530345 0.73462394 1.91645527
#> [2,] -0.55429733 1.5236359 -0.44435298 -0.8293970 0.02137105 1.62014297
#> [3,] -0.65277340 -0.2365804 -0.33365748 1.6144115 -1.38882044 1.18979400
#> [4,] 0.06979050 -0.1798988 -0.01576611 -1.6377538 -0.85286458 0.07352177
#> [5,] -0.58450206 -1.7411024 -0.13801223 -0.4512545 -0.14449068 -0.77910581
#> [6,] -1.26045088 0.2087882 -0.72491869 0.6856348 -0.43163916 -1.62517684
#> [7,] 0.03836372 0.4497529 -0.14266131 -0.6019886 1.36552523 0.14761767
#> [8,] -0.52083678 -0.7626927 0.12068723 -1.6146444 0.79777307 2.41326136
#> [9,] 0.30631735 0.9848562 -0.73193720 -1.2236526 -0.59785470 -0.13138414
#> [10,] -1.67366702 0.8709858 -0.80617294 -0.7559627 1.76655800 0.58235105
#> [,7] [,8] [,9] [,10]
#> [1,] 0.23481727 -0.2019611 -1.0254989 0.97524648
#> [2,] -0.11230652 1.7600720 1.4330532 -0.46696929
#> [3,] 1.80315380 -0.7423216 -0.3679811 -0.83315697
#> [4,] -0.30160066 0.6676371 2.0313025 -0.07749897
#> [5,] 0.04835639 0.8040776 -0.2039855 -0.75413152
#> [6,] 0.56503309 -0.1387350 -0.4091602 -1.09688456
#> [7,] 1.17868940 0.5279960 -0.5626160 -0.09706896
#> [8,] -1.64916614 -0.6481176 1.7608488 -0.69320924
#> [9,] -0.44649730 0.4507879 1.4486604 0.60032266
#> [10,] -1.48612450 0.1245139 -0.9288625 -1.10028291Response derivation
A response can now be added to the dataset by the
simulation_Data function. We have to specifiy the support
of the response, i.e. the explanatory variables that will be used in the
linear model created to compute the response. The support is given by
the supp vector whose entries are either \(0\) or \(1\). The length of the supp
vector must be equal to the number of explanatory variables and if the
\(i\)entry is equal to \(1\), it means that the \(i\)variable will be used to derive the
response value, whereas if the \(i\)entry is equal to \(0\), it means that the \(i\)variable will not be used to derive the
response value (beta<-rep(0,length(supp))). The values
of the coefficients for the explanatory variables that are in the
support of the response are random (either absolute value and sign) and
given by
beta[which(supp==1)]<-runif(sum(supp),minB,maxB)*(rbinom(sum(supp),1,.5)*2-1).
Hence, the user can specify their minimal absolute value with the
minB option and their maximal absolute value with the
maxB option. The stn option is a scaling
factor for the noise added to the response vector
((t(beta)%*%var(X)%*%beta)/stn, with X the
matrix of explanatory variables). The higher the stn value,
the smaller the noise: for instance for a given X dataset,
an stn value \(\alpha\)
times larger will result in a noise exactly \(\sqrt{\alpha}\) times smaller.
set.seed(3141)
supp<-c(1,1,1,0,0,0,0,0,0,0)
minB<-1
maxB<-2
stn<-50
firstdataset=simulation_DATA(X,supp,minB,maxB,stn)
firstdataset
#> $X
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.69029599 -1.7343187 0.38993973 0.7530345 0.73462394 1.91645527
#> [2,] -0.55429733 1.5236359 -0.44435298 -0.8293970 0.02137105 1.62014297
#> [3,] -0.65277340 -0.2365804 -0.33365748 1.6144115 -1.38882044 1.18979400
#> [4,] 0.06979050 -0.1798988 -0.01576611 -1.6377538 -0.85286458 0.07352177
#> [5,] -0.58450206 -1.7411024 -0.13801223 -0.4512545 -0.14449068 -0.77910581
#> [6,] -1.26045088 0.2087882 -0.72491869 0.6856348 -0.43163916 -1.62517684
#> [7,] 0.03836372 0.4497529 -0.14266131 -0.6019886 1.36552523 0.14761767
#> [8,] -0.52083678 -0.7626927 0.12068723 -1.6146444 0.79777307 2.41326136
#> [9,] 0.30631735 0.9848562 -0.73193720 -1.2236526 -0.59785470 -0.13138414
#> [10,] -1.67366702 0.8709858 -0.80617294 -0.7559627 1.76655800 0.58235105
#> [,7] [,8] [,9] [,10]
#> [1,] 0.23481727 -0.2019611 -1.0254989 0.97524648
#> [2,] -0.11230652 1.7600720 1.4330532 -0.46696929
#> [3,] 1.80315380 -0.7423216 -0.3679811 -0.83315697
#> [4,] -0.30160066 0.6676371 2.0313025 -0.07749897
#> [5,] 0.04835639 0.8040776 -0.2039855 -0.75413152
#> [6,] 0.56503309 -0.1387350 -0.4091602 -1.09688456
#> [7,] 1.17868940 0.5279960 -0.5626160 -0.09706896
#> [8,] -1.64916614 -0.6481176 1.7608488 -0.69320924
#> [9,] -0.44649730 0.4507879 1.4486604 0.60032266
#> [10,] -1.48612450 0.1245139 -0.9288625 -1.10028291
#>
#> $Y
#> [1] -4.2132936 3.5039588 0.3332549 -0.4924011 -2.5391834 1.8674007
#> [7] 0.6678607 -0.4589311 0.6353867 3.8091855
#>
#> $support
#> [1] 1 1 1 0 0 0 0 0 0 0
#>
#> $beta
#> [1] -1.754996 1.964992 1.041431 0.000000 0.000000 0.000000 0.000000
#> [8] 0.000000 0.000000 0.000000
#>
#> $stn
#> [1] 50
#>
#> $sigma
#> [,1]
#> [1,] 0.3493447
#>
#> attr(,"class")
#> [1] "simuls"Multiple datasets and checks
To generate multiple datasets, repeat steps 2 and 3, for instance use
a for loop. We create \(NDatasets=200\) datasets and assign them to
the objects DATA_exemple1_nb_1 to
DATA_exemple1_nb_200.
set.seed(3141)
NDatasets=200
for(i in 1:NDatasets){
X<-simulation_X(N,CM)
assign(paste("DATA_exemple1_nb_",i,sep=""),simulation_DATA(X,supp,minB,maxB,stn))
}We now check the correlation structure of the explanatory variable. First we compute the mean correlation matrix.
corr_sum=matrix(0,length(group),length(group))
for(i in 1:NDatasets){
corr_sum=corr_sum+cor(get(paste("DATA_exemple1_nb_",i,sep=""))$X)
}
corr_mean=corr_sum/NDatasetsThen we display and plot that the mean correlation matrix.
corr_mean
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.0000000000 -0.0008611262 0.0193872629 0.0192496952 -0.012147407
#> [2,] -0.0008611262 1.0000000000 -0.0520800766 0.0144798781 0.006237499
#> [3,] 0.0193872629 -0.0520800766 1.0000000000 0.0008693002 -0.021373842
#> [4,] 0.0192496952 0.0144798781 0.0008693002 1.0000000000 0.007753693
#> [5,] -0.0121474071 0.0062374994 -0.0213738420 0.0077536931 1.000000000
#> [6,] -0.0089756967 -0.0404111300 0.0344817040 0.0081889675 0.018018674
#> [7,] -0.0082911544 0.0072612885 -0.0233188445 -0.0380192689 0.023833224
#> [8,] 0.0272233550 -0.0066654749 -0.0487035643 0.0172624295 0.043181249
#> [9,] -0.0145986545 0.0071146338 0.0364868095 -0.0020153080 -0.027733046
#> [10,] 0.9422544272 -0.0071281448 0.0264886880 0.0221950354 -0.003811061
#> [,6] [,7] [,8] [,9] [,10]
#> [1,] -0.008975697 -0.008291154 0.027223355 -0.014598655 0.942254427
#> [2,] -0.040411130 0.007261289 -0.006665475 0.007114634 -0.007128145
#> [3,] 0.034481704 -0.023318845 -0.048703564 0.036486809 0.026488688
#> [4,] 0.008188968 -0.038019269 0.017262430 -0.002015308 0.022195035
#> [5,] 0.018018674 0.023833224 0.043181249 -0.027733046 -0.003811061
#> [6,] 1.000000000 -0.015449494 -0.004054573 0.006159349 0.003444504
#> [7,] -0.015449494 1.000000000 -0.002105066 0.005052182 -0.018230902
#> [8,] -0.004054573 -0.002105066 1.000000000 -0.003169857 0.021688766
#> [9,] 0.006159349 0.005052182 -0.003169857 1.000000000 -0.013388952
#> [10,] 0.003444504 -0.018230902 0.021688766 -0.013388952 1.000000000
plot(abs(corr_mean))
coef_sum=rep(0,length(group))
names(coef_sum)<-paste("x",1:length(group),sep="")
error_counter=0
for(i in 1:NDatasets){
tempdf=data.frame(cbind(Y=get(paste("DATA_exemple1_nb_",i,sep=""))$Y,
get(paste("DATA_exemple1_nb_",i,sep=""))$X))
tempcoef=coef(lm(Y~.-1,data=tempdf))
if(is.null(tempcoef)){
cat("Error in lm fit, skip coefficients\n")
error_counter=error_counter+1
} else{
coef_sum=coef_sum+abs(tempcoef)
}
}
error_counter
#> [1] 0
coef_mean=coef_sum/NDatasetsAll fits were sucessful. Then we display and plot that the mean coefficient vector values.
coef_mean
#> x1 x2 x3 x4 x5 x6
#> 1.508327e+00 1.518967e+00 1.491694e+00 1.437883e-15 1.857978e-15 2.524863e-15
#> x7 x8 x9 x10
#> 1.848380e-15 2.314167e-15 2.601342e-15 6.269818e-15
barplot(coef_mean)
abline(h=(minB+maxB)/2,lwd=2,lty=2,col="blue")
Reduce the noise in the response for the new responses by a factor \(\sqrt{5000/50}=10\). \(1/stn\cdot \beta_{support}^t\mathrm{Var}(X)\beta_{support}\) where \(\beta_{support}\) is the vector of coefficients wh
set.seed(3141)
stn <- 5000
for(i in 1:NDatasets){
X<-simulation_X(N,CM)
assign(paste("DATA_exemple1_bis_nb_",i,sep=""),simulation_DATA(X,supp,minB,maxB,stn))
}Since it is the same explanatory dataset for response generation, we can compare the \(\sigma\) between those \(NDatasets=200\) datasets.
stn_ratios=rep(0,NDatasets)
for(i in 1:NDatasets){
stn_ratios[i]<-get(paste("DATA_exemple1_nb_",i,sep=""))$sigma/
get(paste("DATA_exemple1_bis_nb_",i,sep=""))$sigma
}
all(sapply(stn_ratios,all.equal,10))
#> [1] TRUEAll the ratios are equal to 10 as anticipated.
Since, the correlation structure is the same as before, we do not
need to check it again. As befor, we infer the coefficients values of a
linear model using the lm function.
coef_sum_bis=rep(0,length(group))
names(coef_sum_bis)<-paste("x",1:length(group),sep="")
error_counter_bis=0
for(i in 1:NDatasets){
tempdf=data.frame(cbind(Y=get(paste("DATA_exemple1_bis_nb_",i,sep=""))$Y,
get(paste("DATA_exemple1_bis_nb_",i,sep=""))$X))
tempcoef=coef(lm(Y~.-1,data=tempdf))
if(is.null(tempcoef)){
cat("Error in lm fit, skip coefficients\n")
error_counter_bis=error_counte_bisr+1
} else{
coef_sum_bis=coef_sum_bis+abs(tempcoef)
}
}
error_counter_bis
#> [1] 0
coef_mean_bis=coef_sum_bis/NDatasetsAll fits were sucessful. Then we display and plot that the mean
coefficient vector values. As expected, the noise reduction enhances the
retrieval of the true mean coefficient absolute values by
the models.
coef_mean_bis
#> x1 x2 x3 x4 x5 x6
#> 1.508327e+00 1.518967e+00 1.491694e+00 1.437883e-15 1.857978e-15 2.524863e-15
#> x7 x8 x9 x10
#> 1.848380e-15 2.314167e-15 2.601342e-15 6.269818e-15
barplot(coef_mean_bis)
abline(h=(minB+maxB)/2,lwd=2,lty=2,col="blue")
The simulation process looks sucessfull. What are the confidence indices for those variables?
