Using R with examining apartment house energy consumption effectiveness (DEA effectiveness)

This example is test only. Data itself is “real” data and I think you could do that example also with own data if available. In this example I used library benchmarking.

First set the library and  download data from source

#################################
library(Benchmarking)
#################################
#DEA 1
#output haitakkeet huoneiston suhteen vesi, lämpö ja sähkö
#input vuokransuuruus, asuntojen lukumäärä
#/*lukee suoraan sas-ohjelman csv exportatun datan*/
deamk <- read.table(“http://ekqvist.goeuropeinfo.com/rbloggerqvist/data/deamk.csv”, header=TRUE, sep=”;”, na.strings=”NA”, dec=”,”,
strip.white=TRUE)

#dea1
x1 <- with(deamk, cbind(vuokra_e_nelio_kk, asuntoja)) #INPUT PANOKSET
x1 <- with(deamk, cbind(asuntoja, vuokra_e_nelio_kk)) #INPUT PANOKSET

x1 <- with(deamk, cbind(vuokra_e_nelio_kk)) #INPUT PANOKSET
y1 <- with(deamk, cbind(vesi_vrk_as_litraa, sahko_MWh, lampo_MWh)) # OUTPUT TUOTOKSET

dea.plot.frontier(x1,y1,txt=TRUE, main=”DEA Tehokkuusrintama HeKa aineistosta 2011″)
dea.plot(x1,y1)

#If you want graph into local folder change appropriate folder name below
# example – output graph to jpeg file
jpeg(“i:/todo/asuntodata_dea/deamk.jpg”)
dea.plot.frontier(x1,y1,txt=TRUE, main=”DEA Tehokkuusrintama HeKa aineistosta 2011″)
dev.off()

#note before that install “benchmarking” library
e_dea2011out <- dea(x1,y1, RTS=”vrs”, ORIENTATION=”out”)
e_dea2011in <- dea(x1,y1, RTS=”vrs”, ORIENTATION=”in”)   #TÄTÄ KÄYTETTY ANALYYSISITEKSTISSÄ

#saving results into file
#tehokkuusluvut tiedostoon
write.table(eff(e_dea2011out), file = “i:/todo/asuntodata_dea/e_dea2011out.txt”, sep = “\t”, col.names = NA, qmethod = “double”)
write.table(eff(e_dea2011in), file = “i:/todo/asuntodata_dea/e_dea2011in.txt”, sep = “\t”, col.names = NA, qmethod = “double”)

#..more informative by using summary
sink(“i:/todo/asuntodata_dea/summary_deamk.txt”)
summary(e_dea2011in)
sink()

#plotting effectiveness data with density graph
jpeg(“i:/todo/asuntodata_dea/dea1_density_mkkiinteistot.jpg”)
eff.dens.plot(e_dea2011in , main=”DEA tiheysjakauma HeKa 2011″)
dev.off()
#and same pic here


#070512

#we are also interested which variable explain most  best effectiveness..

#Meitä kiinnostaa myös mitkä valituista muuttujista selittävät parhaiten
# tehokkuutta
# vesi, lämpö ja sähkö, vuokransuuruus
#testataan tätä lineaarisella regressiolla
# yhdistetään eff tiedostoon alkuperäiset muuttujat
#luetaan eff tiedosto sellainsenaan edellä ajetusta

#here we read above effectiveness data table (note 1= best unit)
effvalmis <- read.table(“http://ekqvist.goeuropeinfo.com/rbloggerqvist/data/eff_dea2011valmis.txt”, header=TRUE, sep=”\t”, na.strings=”NA”, dec=”.”, strip.white=TRUE)
#in this operation we could use merge and unique link variable “id”
yhdista1 <- merge(deamk, effvalmis, by.x = “id”, by.y = “id”, all = TRUE)
attach(yhdista1)
names(yhdista1)
summary(yhdista1)

#Examine variables more details – TUTKITAAN MUUTTUJIA TARKEMMIN
hist(yhdista1$x)
hist(yhdista1$vuokra_e_nelio_kk)
hist(yhdista1$vesi_vrk_as_litraa)
hist(yhdista1$sahko_MWh)
hist(yhdista1$lampo_MWh)

#TUTKITAAN MUUTTUJIA TARKEMMIN
jpeg(“j:/todo/asuntodata_dea/hist_matriisikuva.jpg”)
layout(matrix(1:6,2,3))
hist(yhdista1$x)
hist(yhdista1$vuokra_e_nelio_kk)
hist(yhdista1$vesi_vrk_as_litraa)
hist(yhdista1$sahko_MWh)
hist(yhdista1$lampo_MWh)
dev.off()

#DENSITY
jpeg(“j:/todo/asuntodata_dea/density_matriisikuva.jpg”)
layout(matrix(1:6,2,3))
plot(density(yhdista1$x,na.rm=TRUE))
plot(density(yhdista1$vuokra_e_nelio_kk,na.rm=TRUE))
plot(density(yhdista1$vesi_vrk_as_litraa,na.rm=TRUE))
plot(density(yhdista1$sahko_MWh,na.rm=TRUE))
plot(density(yhdista1$lampo_MWh,na.rm=TRUE))
dev.off()

#SORTED VERSIO
jpeg(“j:/todo/asuntodata_dea/sorted_matriisikuva.jpg”)
layout(matrix(1:6,2,3))
plot(sort(yhdista1$x),pch=”1″)
plot(sort(yhdista1$vuokra_e_nelio_kk),pch=”2″)
plot(sort(yhdista1$vesi_vrk_as_litraa),pch=”3″)
plot(sort(yhdista1$sahko_MWh),pch=”4″)
plot(sort(yhdista1$lampo_MWh),pch=”5″)
dev.off()


#valitaan vain tietyt muuttujat pairs kuvaajaan
yhdista2 <- with(yhdista1, cbind(x, vuokra_e_nelio_kk, vesi_vrk_as_litraa, sahko_MWh, lampo_MWh)) #INPUT PANOKSET

#pairs
jpeg(“j:/todo/asuntodata_dea/pairs_matriisikuva.jpg”)
pairs(yhdista2)
dev.off()

#selitetään tehokkuuskerrointa panos ja tuotosmuuttujilla
sink(“i:/todo/asuntodata_dea/summary_deamk_fit_kaikki_muuttujat.txt”)
lm.r = lm(x ~ vuokra_e_nelio_kk + vesi_vrk_as_litraa + sahko_MWh + lampo_MWh)
summary(lm.r)
sink()

 

#vain tuotosmuuttujat
sink(“i:/todo/asuntodata_dea/summary_deamk_fit_vain_tuotos_muuttujat.txt”)
lm0.r = lm(x ~ vesi_vrk_as_litraa + sahko_MWh + lampo_MWh)
summary(lm0.r)
sink()

 

 

#kertoimet
coef(lm.r) # gives the model’s coefficients

#jäännöstermit eli kuinka paljon annetuilla arvoilla ennuste poikkeaa toteumasta
resid(lm.r) # gives the residual errors in Y

fitted(lm.r) # gives the predicted values for Y

#RESIDUAALIEN TUTKIMINEN

layout(matrix(1:4,2,2))
plot(lm.r)

#KAIKKI
jpeg(“i:/todo/asuntodata_dea/dea_residuaalitutkinta_mkkiinteistot_kaikki.jpg”)
layout(matrix(1:4,2,2))
plot(lm.r)
dev.off()

#VAIN TUOTOSMUUTTUJAT
jpeg(“i:/todo/asuntodata_dea/dea_residuaalitutkinta_mkkiinteistot_tuotosmuuttujat.jpg”)
layout(matrix(1:4,2,2))
plot(lm0.r)
dev.off()

#VUOKRA
jpeg(“i:/todo/asuntodata_dea/dea_residuaalitutkinta_mkkiinteistot_vuokra.jpg”)
layout(matrix(1:4,2,2))
plot(lm1.r)
dev.off()

#VESI
jpeg(“i:/todo/asuntodata_dea/dea_residuaalitutkinta_mkkiinteistot_vesi.jpg”)
layout(matrix(1:4,2,2))
plot(lm2.r)
dev.off()

#SÄHKÖ
jpeg(“i:/todo/asuntodata_dea/dea_residuaalitutkinta_mkkiinteistot_sahko.jpg”)
layout(matrix(1:4,2,2))
plot(lm3.r)
dev.off()

#LÄMPÖ
jpeg(“i:/todo/asuntodata_dea/dea_residuaalitutkinta_mkkiinteistot_lampo.jpg”)
layout(matrix(1:4,2,2))
plot(lm4.r)
dev.off()

#below in Finnish only

ennen kuin regressiomallin tulokset voisi hyväksyä, tulisi arvioida sen sopivuus
selittää dataa. yksi tapa toteuttaa tämä on tutkia residuaaleja visuaalisesti.
mikäli malli on soveltuva, tällöin residuaalivirheet tulisi olla satunnaisia ja
normaalisti jakautuneita.

Lisäksi, poistamalla “tapaus” one case ei tulisi merkittävästi vaikuttaa
mallin sopivuuteen. R mahdollistaa neljän graafisen esityksen tämän arvoimiseksi
käyttämällä plot komentoa.

upper left corner /vasen yläkulma
residuaalit vs mallilla sovitetut arvot. residuaalien tulisi olla satunnaisesti jakautuneita keskiviivan ympärillä. keskiviiva kuvaa virhearvoa nolla. pisteiden suhteen ei siis tulisi olla merkittävää hajontaa tämän viivan ympärillä.

lower left corner / vasen alakulma
standardi Q-Q kuvaaja, tämän kuvaajan tulisi vihjata ovatko residuaalivirheet
normaalisti jakautuneita.

upper right corner /oikea yläkulma
scale-location standardisoitujen residuaalien neliöjuuri (virhetermien tai virhee neliöjuuri) sovituksen tai mallin arvojen funktiona. tässäkin tapauksessa tässä ei tulisi olla tai nähtävissä selkeää trendiä. mikäli on, niin malli tulisi tarkistaa tai mahdollisesti poistaa outlier havaintoja.

lower right corner / oikea alakulma
näyttää jokaisen havainnon mitan eli  “vipuvaikutuksen” tai “voiman” kun määritetään regression tuottamia tuloksia. Cooks etäisyys kuvataan kahdella viivalla, joka on eräs mittari kuvaamaan havaintojen merkitsevyyttä regressiomallissa. pieni etäisyys tarkoittaa sitä, että havainnon / havaintojen poistamisella on vain pieni vaikutus regression tuloksiin. Etäisyydet, jotka ylittävät arvon yksi ovat kyseenalaisia ja viittaavat aineistossa mahdollisesti oleviin outliereihin tai muuten vääränlaiseen malliin (malliasetelmaan).

#Using the Results of a Regression to Make Predictions

vuokra_e_nelio_kk_new=c(3,4,4.5,4.6,4.5,5.5,6.9,5.8,3.6,4.7,5.3,6.2,5.1,3.9,4.3);vuokra_e_nelio_kk_new
predict(lm.r,data.frame(vuokra_e_nelio_kk = vuokra_e_nelio_kk_new), level = 0.9, interval = “confidence”)
where ‘lwr’ is the lower limit of the confidence interval and ‘upr’ is the upper limit of the confidence interval. R does not contain a feature for finding the confidence intervals for predicted values of the independent variable for specified values of dependent variables.

# data plotting as scatter plot/ pistekuvaaja näitä voidaan tehdä niin monta kuin on muuttujia

#SCATTER PLOT

jpeg(“j:/todo/asuntodata_dea/scatterplot_withoutfitline_rent.jpg”)
plot(x, vuokra_e_nelio_kk) dev.off()
plot(x, vesi_vrk_as_litraa)
plot(x, sahko_MWh)
plot(x, lampo_MWh)
plot(x, asuntoja)

 

#FITLINE / KUVAAJA FIT JA TOTEUTUNEET
lm1.r = lm(x ~ vuokra_e_nelio_kk)
plot(x ~ vuokra_e_nelio_kk, yhdista1)
abline(lm1.r$coef,lty=5)

lm2.r = lm(x ~ vesi_vrk_as_litraa)
plot(x ~ vesi_vrk_as_litraa, yhdista1)
abline(lm2.r$coef,lty=5)

lm3.r = lm(x ~ sahko_MWh)
plot(x ~ sahko_MWh, yhdista1)
abline(lm3.r$coef,lty=5)

lm4.r = lm(x ~ lampo_MWh )
plot(x ~ lampo_MWh, yhdista1)
abline(lm4.r$coef,lty=5)

lm5.r = lm(x ~ asuntoja)
plot(x ~ asuntoja, yhdista1)
abline(lm5.r$coef,lty=5)

# example – output graph to jpeg file
jpeg(“i:/todo/asuntodata_dea/fitline_vuokra.jpg”)
lm1.r = lm(x ~ vuokra_e_nelio_kk)
plot(x ~ vuokra_e_nelio_kk, yhdista1)
abline(lm1.r$coef,lty=5)
dev.off()

# example – output graph to jpeg file
jpeg(“i:/todo/asuntodata_dea/fitline_vesi.jpg”)
lm2.r = lm(x ~ vesi_vrk_as_litraa)
plot(x ~ vesi_vrk_as_litraa, yhdista1)
abline(lm2.r$coef,lty=5)
dev.off()

# example – output graph to jpeg file
jpeg(“i:/todo/asuntodata_dea/fitline_sahko.jpg”)
lm3.r = lm(x ~ sahko_MWh)
plot(x ~ sahko_MWh, yhdista1)
abline(lm3.r$coef,lty=5)
dev.off()

# example – output graph to jpeg file
jpeg(“i:/todo/asuntodata_dea/fitline_lampo.jpg”)
lm4.r = lm(x ~ lampo_MWh )
plot(x ~ lampo_MWh, yhdista1)
abline(lm4.r$coef,lty=5)
dev.off()

# example – output graph to jpeg file
jpeg(“i:/todo/asuntodata_dea/fitline_asuntoja.jpg”)
lm5.r = lm(x ~ asuntoja)
plot(x ~ asuntoja, yhdista1)
abline(lm5.r$coef,lty=6)
dev.off()

#IN THE SAME PICTURE / SAMAAN KUVAAN
jpeg(“i:/todo/asuntodata_dea/fitline_matriisikuva.jpg”)
layout(matrix(1:6,2,3))
lm1.r = lm(x ~ vuokra_e_nelio_kk)
plot(x ~ vuokra_e_nelio_kk, yhdista1)
abline(lm1.r$coef,lty=5)
lm2.r = lm(x ~ vesi_vrk_as_litraa)
plot(x ~ vesi_vrk_as_litraa, yhdista1)
abline(lm2.r$coef,lty=5)
lm3.r = lm(x ~ sahko_MWh)
plot(x ~ sahko_MWh, yhdista1)
abline(lm3.r$coef,lty=5)
lm4.r = lm(x ~ lampo_MWh )
plot(x ~ lampo_MWh, yhdista1)
abline(lm4.r$coef,lty=5)
lm5.r = lm(x ~ asuntoja)
plot(x ~ asuntoja, yhdista1)
abline(lm5.r$coef,lty=5)
dev.off()

#MUUTTUJAT SUHTEESSA TEHOKKUUSKERTOIMEEN
#LISÄKSI PIIRRETÄÄN KUNKIN MUUTTUJAN FIT-VIIVA
jpeg(“j:/todo/asuntodata_dea/fitline_kaikki_tuotosmuuttujat_samassakuvassa2.jpg”)
matplot(cbind(lampo_MWh,sahko_MWh,vesi_vrk_as_litraa),x,xlab=”x”,ylab=”y”)
#matplot(cbind(x,(lampo_MWh,sahko_MWh,vesi_vrk_as_litraa),xlab=”x”,ylab=”y”)
abline(lm4.r,lty=1)#LÄMPÖ
abline(lm3.r,lty=2)#SÄHKÖ
abline(lm2.r,lty=3)#VESI
dev.off()

Figure : Original x shown with “1”, with small error as “2” and with large error as “3”. The regression lines for the no measurement error, small error and large error are
shown as solid, dotted and dashed lines respectively.

#DENSITY / JAKAUMAT
jpeg(“i:/todo/asuntodata_dea/density_matriisikuva.jpg”)
layout(matrix(1:4,2,2))
plot(density(lampo_MWh,na.rm=TRUE))
plot(density(sahko_MWh,na.rm=TRUE))
plot(density(vesi_vrk_as_litraa,na.rm=TRUE))
dev.off()

#EXAMINING CORRELATION / KATSOTAAN VIELÄ KORRELAATIOT
data <- with(yhdista1, cbind(x, asuntoja, vuokra_e_nelio_kk, vesi_vrk_as_litraa, sahko_MWh, lampo_MWh))

# First Correlogram Example
jpeg(“i:/todo/asuntodata_dea/correlogram1_.jpg”)
library(corrgram)
corrgram(data, order=TRUE, lower.panel=panel.shade,
upper.panel=panel.pie, text.panel=panel.txt,
main=”DEA tehokkuuskertoimen (x) ja mallin muuttujien välinen korrelaatio”)
dev.off()

# Second Correlogram Example
# library(corrgram)
jpeg(“i:/todo/asuntodata_dea/correlogram2_.jpg”)
corrgram(data, order=TRUE, lower.panel=panel.ellipse,
upper.panel=panel.pts, text.panel=panel.txt,
diag.panel=panel.minmax,
main=”DEA tehokkuuskertoimen (x) ja mallin muuttujien välinen korrelaatio”)
dev.off()

# Third Correlogram Example
jpeg(“i:/todo/asuntodata_dea/correlogram3_.jpg”)
library(corrgram)
corrgram(data, order=NULL, lower.panel=panel.shade,
upper.panel=NULL, text.panel=panel.txt,
main=”DEA tehokkuuskertoimen (x) ja mallin muuttujien välinen korrelaatio”)
dev.off()

library(psych)
library(corrgram)

#PIIRRETÄÄN KUVAA
jpeg(“i:/todo/asuntodata_dea/matrix_scatterplots_.jpg”)
pairs(data, main=”DEA tehokkuuskertoimen (x) ja mallin muuttujien välinen korrelaatio”)   #draw a matrix of scatter plots for the first 5 variables
dev.off()

jpeg(“i:/todo/asuntodata_dea/histograms_for_all_variables_.jpg”)
multi.hist(data) #draw the histograms for all variables
dev.off()

jpeg(“i:/todo/asuntodata_dea/scatter_histograms_correlations_for_all_variables_.jpg”)
pairs.panels(data) #draws scatterplots, histograms, and shows correlations  (source found in the psych package)
dev.off()

That’s it…..

Copyright MySci 2018
Tech Nerd theme designed by Siteturner