diff --git a/Census2022-CER-mixed_model_0910.R b/Census2022-CER-mixed_model_0910.R
new file mode 100644
index 0000000000000000000000000000000000000000..871534dd163310a5b7ff3bc929e2ef5b302af6ee
--- /dev/null
+++ b/Census2022-CER-mixed_model_0910.R
@@ -0,0 +1,435 @@
+####################################################################################
+## Data Analysis for ESRC Transformative project
+## - Using the Commission for Energy Regulation (CER)'s Irish Smart Meter Trial data
+##   - http://www.ucd.ie/issda/data/commissionforenergyregulationcer/
+
+## This file is used to produce midweek (Tues-Thurs) and weekend (Sat-Sun) mixed models which test the ability of 
+## various power (electricity consumption) variables to predict the household attributes of interest 
+
+## The energy data is for 28 days (4 weeks) of Oct 2009 & 2010
+## October 2009 was before the smart meter trials started and
+## was also close to the pre-trial survey date
+
+## Comparativ emodels also run for 2010
+
+## For the latest version of this code go to: https://github.com/dataknut/Census2022
+
+## This work was funded by RCUK through the ESRC's Transformative Social Science Programme via the
+## "Census 2022: Transforming Small Area Socio-Economic Indicators through 'Big Data'" Project 
+## - http://gtr.rcuk.ac.uk/project/2D2CD798-4F04-4399-B1AF-D810A233DD21
+## - http://www.energy.soton.ac.uk/tag/census2022/
+ 
+## Copyright (C) 2014  University of Southampton
+
+## Author: Sharon Lin (X.Lin@soton.ac.uk) 
+##  [Energy & Climate Change, Faculty of Engineering & Environment, University of Southampton]
+
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License 
+## (http://choosealicense.com/licenses/gpl-2.0/), or (at your option) any later version.
+
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+
+## #YMMV - http://en.wiktionary.org/wiki/YMMV
+
+
+rm(list=ls())
+library(lme4)
+
+
+
+## load in energy data
+setwd("//soton.ac.uk/ude/PersonalFiles/Users/xl25g12/mydocuments/census2022/CER_analysis/data")   ## the file is in working order
+## egdata = read.table("OctHH_midwk_long", header=T, stringsAsFactors=FALSE)   ## read in R exported txt file, no need sep="t"; long form of mid week data for all Oct 2009 and 2010 data
+egdata = read.table("CER_OctHH_wkend_long", header=T, stringsAsFactors=FALSE)   ## read in R exported txt file, no need sep="t"; long form of mid week data for all Oct 2009 and 2010 data
+eg09 = egdata[egdata$DateOct<365,] 
+eg10 = egdata[egdata$DateOct>365,] 
+
+## load in HH survey data, note survey content for pre-trial and post-trial are different
+data1 = read.table("HH09_pretrial.txt",header=T,stringsAsFactors=FALSE,sep="\t")   ## read in txt file, need sep="t", HHpresurvey data
+HH09FA = read.table("HH_floor_pretrial_survey.txt",header=T,stringsAsFactors=FALSE,sep="\t")   ## read in txt file, need sep="t"; pretrial floor area data
+## HH10 = read.table("HH10_posttrial.txt",header=T,stringsAsFactors=FALSE,sep="\t")   ## read in txt file, need sep="t", HHpresurvey data
+HH09 = merge(data1,HH09FA)
+
+## merge energy and HH data together
+## data for 09 and 10 separately
+mydata09=merge(eg09,HH09)
+##mydata10=merge(eg10,HH10)
+
+
+
+
+
+########## mixed effects models are run on these separate files  ###################
+## check on mixed model with income variables using 09 data
+## make a comparison with 0910 data results using 0910
+
+mydata09$dummy_HHchild = NA
+mydata09 = within(mydata09,dummy_HHchild[HHtype_y09==1]<-0)   
+mydata09 = within(mydata09,dummy_HHchild[HHtype_y09==2]<-0)   
+mydata09 = within(mydata09,dummy_HHchild[HHtype_y09==3]<-1)   
+ddply(mydata09, .(dummy_HHchild),nrow)
+ddply(mydata09, .(HHtype_y09),nrow)
+
+## make a new dummy vriable column D_unemployed 
+mydata09$dummy_unemp_R = NA
+mydata09 = within(mydata09,dummy_unemp_R[employment_y09==1]<-0)   
+mydata09 = within(mydata09,dummy_unemp_R[employment_y09==4]<-1)   
+mydata09 = within(mydata09,dummy_unemp_R[employment_y09==5]<-1)   
+mydata09 = within(mydata09,dummy_unemp_R[employment_y09==2]<-NA)   
+mydata09 = within(mydata09,dummy_unemp_R[employment_y09==3]<-NA)   
+mydata09 = within(mydata09,dummy_unemp_R[employment_y09==6]<-NA)
+
+ddply(mydata09, .(dummy_unemp_R),nrow)
+
+## make a second new dummy vriable column dum_em_mixed,  (1-3)
+
+
+
+## prepare dummy variable for "number of residents"
+mydata09["dum_noResident"] = NA
+mydata09 = within(mydata09,dum_noResident[HHtype_y09==1]<-0)   
+mydata09 = within(mydata09,dum_noResident[HHtype_y09==3&HHchild_y09==1]<-0)   
+mydata09 = within(mydata09,dum_noResident[HHtype_y09==2&HHadult_y09==1]<-0)   
+mydata09 = within(mydata09,dum_noResident[HHadult_y09>1]<-1)   
+mydata09 = within(mydata09,dum_noResident[HHchild_y09>1]<-1)   
+ddply(mydata09, .(dum_noResident),nrow)
+names(mydata09)
+
+
+## change income to approximated mid ranged numbers
+mydata09$income_num = mydata09$income   ## add a new column
+mydata09 = within(mydata09,income_num[income_num==1]<-7500)   
+mydata09 = within(mydata09,income_num[income_num==2]<-22500)   
+mydata09 = within(mydata09,income_num[income_num==3]<-40000)
+mydata09 = within(mydata09,income_num[income_num==4]<-62500)
+mydata09 = within(mydata09,income_num[income_num==5]<-92500)   
+mydata09 = within(mydata09,income_num[income_num==6]<-NA)   
+
+## change dataframe structure
+mydata09 = within(mydata09,HHtype_y09<-as.factor(HHtype_y09))
+mydata09 = within(mydata09,dummy_unemp_R<-as.factor(dummy_unemp_R))   
+mydata09 = within(mydata09,dummy_HHchild<-as.factor(dummy_HHchild))
+mydata = mydata09
+
+
+## the following code are used for weekend data 
+final =0
+for (i in c(11,9,12, 3, 5, 7,16,17) )
+{
+  energy  = mydata[,i]
+  
+  ## C1:  test = lmer(energy  ~ 1 + (1|HHID) + num_resident + log(income_num) + emp_dummy, data=mydata, na.action="na.omit")  ## RE on intercept only   
+  ## C2:  test = lmer(energy  ~ 1 + (1|HHID) + num_resident + log(income_num), data=mydata, na.action="na.omit")  ## RE on intercept only 
+  ## C3: test = lmer( energy  ~ 1 + (1|HHID) + num_resident + log(income_num) + HHtype , data=mydata, na.action="na.omit")  ## RE on intercept only 
+  ## C4: test = lmer( energy  ~ 1 + (1|HHID) + num_resident + log(income_num) + dummy_HHchild , data=mydata, na.action="na.omit")  ## RE on intercept only 
+  ## C5: test = lmer(energy ~ 1 + (1|HHID) + num_resident + log(income_num) + dummy_HHchild + dummy_unemp_R , data=mydata, na.action="na.omit")  ## RE on intercept only 
+  ## C4: test = lmer(energy  ~ 1 + (1|ID) + dum_noResident + log(income_num) + dummy_HHchild, data=mydata, na.action="na.omit")  ## RE on intercept only   
+  ## C3: test = lmer(energy  ~ 1 + (1|ID) + dum_noResident + log(income_num), data=mydata, na.action="na.omit")  ## RE on intercept only   
+  ## C1: 
+  test = lmer(energy  ~ 1 + (1|ID) + dum_noResident + log(income_num)+dummy_unemp_R, data=mydata, na.action="na.omit")  ## RE on intercept only   
+   
+  ## Extract fixed effect coefficients
+  out_FE = as.data.frame(cbind(as.data.frame(coef(summary(test)))[,1],as.data.frame(coef(summary(test)))[,3]))
+  colnames(out_FE)[1] = c(get("i"))
+  colnames(out_FE)[2] = c("zsta")
+  
+  ## for 1 RE only: Extract random effect varaince and percentage of ICC (intercorrelation between RE in intercept and residuals)
+    RE1 = as.data.frame(VarCorr(test))$vcov[1]    ## variance of HHID ##  as.data.frame(VarCorr(test))$sdcor[1]  for std. dev.
+    RE2 = as.data.frame(VarCorr(test))$vcov[2]    ## varaince of residual ## 
+  # alternative code
+  ##  as.numeric(summary(test)@REmat[1,3])    ## variance of HHID;  
+  ##  as.numeric(summary(test)@REmat[1,4])    ## std error of HHID
+  ## as.numeric(summary(test)@REmat[2,4])    ## std error of residuals
+  ## k = 2
+  ## tmp11 = RE1 / (RE1 + RE2)    ## % of RE effects in HH
+  ## RE3 = as.numeric(format(round(tmp11,k),nsmall=k))     ## % of RE effects in HH      
+  ## RE = NA   ## REtmp and out are of different length
+  ## out = cbind(out_FE,RE)
+  ## out$RE[1]=RE1
+  ## out$RE[2]=RE2
+  ##  out$RE[3]=RE3
+  
+  testlmer = as.vector(test@pp$X %*% fixef(test))   ## fitted values    ## length of fitted values is 4099, no adjustment is made
+   VarF=var(testlmer)
+  ## VarV = variance (intercept) + varaince(Slope) + 2 cov.correlation * stdev(intercept)stdev(slope)   ## variation generated from RE
+  ##  RE_sdcor =  as.data.frame(VarCorr(test))$sdcor  ## stdev and correlation of two RE effects, and stdev of residuals
+  ## VarV = RE_sdcor[1]^2 
+  ##  VarR = RE_sdcor[2]^2      ## variation from residuals
+    VarV=RE1
+    VarR = RE2
+    R2.m = VarF/(VarF + VarV + VarR)   ## Var(FE)/[Var(FE)+Var(RE) + Var(Res)]
+    R2.c = (VarF+VarV)/(VarF + VarV + VarR)
+    R2.r = VarR / (VarF + VarV + VarR) 
+    R2 = NA   ## REtmp and out are of different length
+    out = cbind(out_FE,R2)
+    out$R2[1]=R2.m
+    out$R2[2]=R2.c
+    out$R2[3]=R2.r
+  
+  
+  
+#   ## for 2 REs in the model: extract RE effects from the model
+#   RE_sdcor =  as.data.frame(VarCorr(test))$sdcor  ## stdev and correlation of two RE effects, and stdev of residuals
+#   ##  final_FE = cbind(final_FE,out_FE)
+#   ## ## final_RE = cbind(final_RE,RE_sdcor)
+#   ## R squared for D1-D5 models,  for mixed with RE in intercept and slope
+#   testlmer = as.vector(test@pp$X %*% fixef(test))   ## fitted values    ## length of fitted values is 4099, no adjustment is made
+#   VarF=var(testlmer)
+#   ##  VarV = variance (intercept) + varaince(Slope) + 2 cov.correlation * stdev(intercept)stdev(slope)   ## variation generated from RE
+#   ##  or VarV = RE_sdcor[1]^2 + 2*RE_sdcor[1]*RE_sdcor[2]*RE_sdcor[3] + RE_sdcor[2]^2   ## variation generated from RE
+#   VarR = RE_sdcor[4]^2      ## variation from residuals
+#   R2.m = VarF/(VarF + VarV + VarR)   ## Var(FE)/[Var(FE)+Var(RE) + Var(Res)]
+#   R2.c = (VarF+VarV)/(VarF + VarV + VarR)
+#   R2.r = VarR / (VarF + VarV + VarR) 
+#   R2 = NA   ## REtmp and out are of different length
+#   out = cbind(out_FE,R2)
+#   out$R2[1]=R2.m
+#   out$R2[2]=R2.c
+#   out$R2[3]=R2.r
+#   
+  final = cbind(final,out) 
+}
+
+
+## RE1 working codes are as above.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+###### models with T1 - T2 (without income variables) 
+
+final=0
+for (i in c(3,7,9,11, 12, 16, 17) )    ## mean, q975, ammean, pktime, base, LF  (ECF has been consistent not converging, taken out)
+  ## i = 16, ECF does not converge
+{
+  energy  = (mydata09[,i])
+  ## D1: test = lmer(energy ~ num_resident + income_num.c + dummy_unemp_R + (1+income_num.c|HHID) , data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  ## D2 test = lmer(energy ~ num_resident + income_num.c + (1+income_num.c|HHID), data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  ## D3: test = lmer(energy ~ num_resident + income_num.c + HHtype + (1+income_num.c|HHID) , data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  ##  D4: test = lmer(energy ~ num_resident + income_num.c + dummy_HHchild + (1+income_num.c|HHID) , data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  
+   ## T1 =  lmer(energy  ~ 1 + dum_noResident + dum_em13_46 + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  
+   T2 = lmer(energy  ~ 1 + dum_noResident + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+   ## T2a = lmer( energy  ~ 1 + dum_noResident + dum_em + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+   ## T3 = lmer(energy  ~ 1 + dum_noResident + dum_HHtype2 + dum_HHtype3 + (1|ID),  data=mydata09, na.action="na.omit") 
+  
+  ## T4 = lmer( energy  ~ 1 + dum_noResident + dum_child + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+  ## T4a = lmer( energy  ~ 1 + dum_child + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+  ## T5 = lmer( energy  ~ 1 + dum_noResident + dum_child + dum_em + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+  ## T6 = lmer( energy  ~ 1 + dum_child + dum_noResident + dum_child*dum_noResident + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+  
+  ## T7 = lmer( energy  ~ 1 + dum_child + dum_noResident + dum_child*dum_noResident + dum_em + (1|ID),  data=mydata09, na.action="na.omit")  ## RE on intercept only 
+  
+  T2 = T2
+  
+  ## Extract fixed effect coefficients
+  out = as.data.frame(cbind(as.data.frame(coef(summary(T2)))[,1],as.data.frame(coef(summary(T2)))[,3]))
+  colnames(out)[1] = c(get("i"))
+  colnames(out)[2] = c("zsta")
+  
+  ## for 1 RE only: Extract random effect varaince and percentage of ICC (intercorrelation between RE in intercept and residuals)
+  ##   command:  as.data.frame(VarCorr(test)) gives  variance-covariance of the model: as.data.frame(VarCorr(test))$vcov;  std.error of the model: as.data.frame(VarCorr(test))$sdcor
+  RE1 = as.data.frame(VarCorr(T2))$vcov[1]    ## variance of RE in the intercept "ID" 
+  RE2 = as.data.frame(VarCorr(T2))$vcov[2]    ## varaince of residual ## 
+  
+  testlmer = as.vector(T2@pp$X %*% fixef(T2))   ## fitted values    ## length of fitted values is 4099, no adjustment is made
+  VarF=var(testlmer)
+  ## VarV = variance (intercept) + varaince(Slope) + 2 cov.correlation * stdev(intercept)stdev(slope)   ## variation generated from RE
+  VarV =RE1
+  VarR = RE2      ## variation from residuals
+  R2.m = VarF/(VarF + VarV + VarR)   ## Var(FE)/[Var(FE)+Var(RE) + Var(Res)]
+  R2.c = (VarF+VarV)/(VarF + VarV + VarR)
+  R2.r = VarR / (VarF + VarV + VarR) 
+  out["R2"]=NA
+  out$R2[1] = R2.m
+  out$R2[2] = R2.c
+  out$R2[3] = R2.r
+  final= cbind(final,out)
+}
+
+
+
+
+
+
+
+
+########  end of mixed effects models on sepearte 09 and 10 periods ###################
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#################  beginning of mixed models on both Oct 09 and Oct 10 data  #######################
+
+######## data for 09 and 10 together  #############
+HH0910 = read.table("HH_both_survey", header=T, stringsAsFactors=FALSE)    ## HH info.
+tmpHH09=HH0910[HH0910$index==0,]     ## HH info in Oct09
+tmpHH10=HH0910[HH0910$index==1,]     ## HH info in Oct10
+
+tdata09=merge(eg09,tmpHH09)   ## HH and enregy infoin Oct 09 where HH has both info (HH participated in both pre and post-survey)
+tdata10=merge(eg10,tmpHH10)   ## HH and enregy infoin Oct 10 where HH has both info (HH participated in both pre and post-survey)
+mydata0910 = rbind(tdata09,tdata10)
+
+mydata0910["dum_em"]=NA    ## in survey data, 1 is employed; 4 is unemployed
+mydata0910 = within(mydata0910,dum_em[employment==1]<-0)  
+mydata0910 = within(mydata0910,dum_em[employment==4]<-1)  
+ddply(mydata0910, .(dum_em),nrow)
+names(mydata0910)
+
+mydata0910["dum_child"] = NA
+mydata0910 = within(mydata0910,dum_child[HHtype==1]<-0)   
+mydata0910 = within(mydata0910,dum_child[HHtype==2]<-0)   
+mydata0910 = within(mydata0910,dum_child[HHtype==3]<-1)   
+ddply(mydata0910, .(dum_child),nrow)
+ddply(mydata0910, .(HHtype),nrow)
+names(mydata0910)
+
+## prepare the new columne,  count for the number of resdients (0/1), 
+mydata0910["dum_noResident"] = NA
+mydata0910 = within(mydata0910,dum_noResident[HHtype==1]<-0)   
+mydata0910 = within(mydata0910,dum_noResident[HHtype==3&HHchild==1]<-0)   
+mydata0910 = within(mydata0910,dum_noResident[HHtype==2&HHadult==1]<-0)   
+mydata0910 = within(mydata0910,dum_noResident[HHadult>1]<-1)   
+mydata0910 = within(mydata0910,dum_noResident[HHchild>1]<-1)   
+ddply(mydata0910, .(dum_noResident),nrow)
+names(mydata0910)
+
+
+##### prepare for dummy_em, with 1-3 value of employment variable as 0, 4-6 as 1  ######
+mydata0910["dum_em13_46"] = NA
+mydata0910 = within(mydata0910,dum_em13_46[employment<4]<-0) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+mydata0910 = within(mydata0910,dum_em13_46[employment>3]<-1)  ## 1 for unemployed people: , 4: unemployed actively seeking work; 5 unemployed not seeking work, 6: retired, 7: carer: looking after relative family)
+ddply(mydata0910, .(dum_em13_46),nrow)
+
+
+## prepare dummy for HHtype, HHtype2 = 1 if HHtype = 2, otherwise 0; HHtype3 = 1 if HHtype = 3, otherwise 0; 
+mydata0910["dum_HHtype2"] = NA
+mydata0910 = within(mydata0910,dum_HHtype2[HHtype==2]<-1) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+mydata0910 = within(mydata0910,dum_HHtype2[HHtype==1]<-0) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+mydata0910 = within(mydata0910,dum_HHtype2[HHtype==3]<-0) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+ddply(mydata0910, .(dum_HHtype2),nrow)
+
+
+mydata0910["dum_HHtype3"] = NA
+mydata0910 = within(mydata0910,dum_HHtype3[HHtype==3]<-1) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+mydata0910 = within(mydata0910,dum_HHtype3[HHtype==1]<-0) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+mydata0910 = within(mydata0910,dum_HHtype3[HHtype==2]<-0) ## 0 for employed people: 1: employee, 2:sel-employed with employees, 3:self-employed without employees;   
+ddply(mydata0910, .(dum_HHtype3),nrow)
+
+
+## change dataframe structure
+mydata0910 = within(mydata0910,MF<-as.factor(MF))   
+mydata0910 = within(mydata0910,HHtype<-as.factor(HHtype))
+mydata0910 = within(mydata0910,age<-as.factor(age))
+mydata0910 = within(mydata0910,employment<-as.factor(employment))   
+mydata0910 = within(mydata0910,socialclass<-as.factor(socialclass))
+mydata0910 = within(mydata0910,HHadult<-as.factor(HHadult))
+mydata0910 = within(mydata0910,dayault<-as.factor(dayault))
+mydata0910 = within(mydata0910,dum_child<-as.factor(dum_child))
+mydata0910 = within(mydata0910,daychild<-as.factor(daychild))
+## mydata0910 = within(mydata0910,bedroom<-as.factor(bedroom))
+## mydata0910 = within(mydata0910,heating<-as.factor(heating))
+## mydata0910 = within(mydata0910,income<-as.factor(income))
+mydata0910 = within(mydata0910,dum_em<-as.factor(dum_em))
+mydata0910 = within(mydata0910,dum_HHtype2<-as.factor(dum_HHtype2))
+mydata0910 = within(mydata0910,dum_HHtype3<-as.factor(dum_HHtype3))
+mydata0910 = within(mydata0910,ID<-as.factor(ID))
+mydata0910 = within(mydata0910,dum_child<-as.factor(dum_child))
+mydata0910 = within(mydata0910,dum_noResident<-as.factor(dum_noResident))
+mydata0910 = within(mydata0910,dum_em13_46<-as.factor(dum_em13_46))
+
+
+########## mixed models for RE effect in the intercept   ############################
+mydata = mydata0910[order(mydata0910$ID),]   ## sort data accoring to HH ID.
+A1 = lmer( dailyammean  ~ 1 + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+A2 = lmer( dailyammean  ~ 1 + dum_noResident + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+A3 = lmer( dailyammean  ~ 1 + dum_noResident + dum_em + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+A4 = lmer( dailyammean  ~ 1 + dum_child + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+A5 = lmer( dailyammean  ~ 1 + dum_child + dum_em + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+A6 = lmer( dailyammean  ~ 1 + dum_child + dum_noResident + dum_child*dum_noResident + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+A7 = lmer( dailyammean  ~ 1 + dum_child + dum_noResident + dum_child*dum_noResident + dum_em + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+
+
+###############  automation to extract tables from RE  ########################## 
+########## loop through one energy variable at a time ##########
+final=0
+for (i in c(3,7,9,11, 12, 17) )    ## mean, q975, ammean, pktime, base, LF  (ECF has been consistent not converging, taken out)
+  ## i = 16, ECF does not converge
+  {
+  energy  = (mydata0910[,i])
+  ## D1: test = lmer(energy ~ num_resident + income_num.c + dummy_unemp_R + (1+income_num.c|HHID) , data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  ## D2 test = lmer(energy ~ num_resident + income_num.c + (1+income_num.c|HHID), data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  ## D3: test = lmer(energy ~ num_resident + income_num.c + HHtype + (1+income_num.c|HHID) , data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  ##  D4: test = lmer(energy ~ num_resident + income_num.c + dummy_HHchild + (1+income_num.c|HHID) , data=mydata, na.action="na.omit")  ## RE on intercept and slope   
+  
+  ## T1 =  lmer(energy  ~ 1 + dum_noResident + dum_em13_46 + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  
+  ## T2 = lmer(energy  ~ 1 + dum_noResident + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  ## T2a = lmer( energy  ~ 1 + dum_noResident + dum_em + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  ## T3 = lmer(energy  ~ 1 + dum_noResident + dum_HHtype2 + dum_HHtype3 + (1|ID),  data=mydata0910, na.action="na.omit") 
+  
+  T4 = lmer( energy  ~ 1 + dum_noResident + dum_child + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  ## T4a = lmer( energy  ~ 1 + dum_child + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  T5 = lmer( energy  ~ 1 + dum_noResident + dum_child + dum_em + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  ## T6 = lmer( energy  ~ 1 + dum_child + dum_noResident + dum_child*dum_noResident + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  
+  ## T7 = lmer( energy  ~ 1 + dum_child + dum_noResident + dum_child*dum_noResident + dum_em + (1|ID),  data=mydata0910, na.action="na.omit")  ## RE on intercept only 
+  
+  T2 = T5
+  
+  ## Extract fixed effect coefficients
+  out = as.data.frame(cbind(as.data.frame(coef(summary(T2)))[,1],as.data.frame(coef(summary(T2)))[,3]))
+  colnames(out)[1] = c(get("i"))
+  colnames(out)[2] = c("zsta")
+  
+  ## for 1 RE only: Extract random effect varaince and percentage of ICC (intercorrelation between RE in intercept and residuals)
+  ##   command:  as.data.frame(VarCorr(test)) gives  variance-covariance of the model: as.data.frame(VarCorr(test))$vcov;  std.error of the model: as.data.frame(VarCorr(test))$sdcor
+  RE1 = as.data.frame(VarCorr(T2))$vcov[1]    ## variance of RE in the intercept "ID" 
+  RE2 = as.data.frame(VarCorr(T2))$vcov[2]    ## varaince of residual ## 
+      
+  testlmer = as.vector(T2@pp$X %*% fixef(T2))   ## fitted values    ## length of fitted values is 4099, no adjustment is made
+  VarF=var(testlmer)
+  ## VarV = variance (intercept) + varaince(Slope) + 2 cov.correlation * stdev(intercept)stdev(slope)   ## variation generated from RE
+  VarV =RE1
+  VarR = RE2      ## variation from residuals
+  R2.m = VarF/(VarF + VarV + VarR)   ## Var(FE)/[Var(FE)+Var(RE) + Var(Res)]
+  R2.c = (VarF+VarV)/(VarF + VarV + VarR)
+  R2.r = VarR / (VarF + VarV + VarR) 
+  out["R2"]=NA
+  out$R2[1] = R2.m
+  out$R2[2] = R2.c
+  out$R2[3] = R2.r
+  final= cbind(final,out)
+}
+
+