I am new in R and it’s my first time using it so I’ll appreciate the help. I am estimating income elasticity for electricity consumption using budget shares. I have data for 8 regions categorized into 5 classes depending on the household size (from A to E). I am regressing budget shares (wi) on log of income(lxi), classes dummy of household size (D) and the interaction between log of income and household class size (lxi*D):

Wi = C + lxi + D + lxi*D + u

I have an unbalanced panel data for 2067 observations saved in .csv format. I attached my data, transformed the Date from factor, made a new data frame including the new date, and finally set the data as a panel data as the code below:

mydata<-read.csv("C:/Users/Fadhila/Desktop/Remeasuring 2015/DataClass-unbalanced.csv", header=T)


date<- as.Date(factor(Date),format= "%m/%d/%Y")


ndata<-plm.data(ndata, index=c("Class", "date"))

I have regressed my model before using Pooled OLS, Fixed (twoways), and Random and both Hausman test and F-test suggested using Random which I doubt so I thought to test for heteroskedasticity and outliers. So I plotted the below model , plotted the leverage versus the residual and compared them with cooks distance:

r <- lm (wi~lxi + Class:lxi, data=ndata)



plot (r)


with(ndata, plot(lxi, cooks.distance(r)))

identify(ndata$lxi, cooks.distance(r))

residuals plot

residual versus cooks distance

Than to estimate how many points are far from the leverage points, but it seems that am doing something wrong as I got all the points to be twice greater than the leverage.

lev = hatvalues(r)




However, I plotted the below to see

plot(ndata$lxi, rstandard(r))


plot(ndata$lxi, lev)

identify(ndata$lxi, lev)


plot(ndata$lxi, ndata$Class)

identify(ndata$lxi, ndata$Class)

enter image description here

enter image description here enter image description here

To double check and see which points influence the results I made the below:

outs <- influencePlot(r)

n <-2

Cooksdist <- as.numeric(tail(row.names(outs[order(outs$CookD), ]), n))

Levr <- as.numeric(tail(row.names(outs[order(outs$Hat), ]), n))

StdRes <- as.numeric(tail(row.names(outs[order(outs$StudRes), ]), n))

plot(ndata$lxi, ndata$wi)

abline(r, col = "blue")

Warning message: In abline(r, col = "blue") : only using the first two of 41 regression coefficients*

points(ndata$lxi[Cooksdist], ndata$wi[Cooksdist], col = "red", pch = 0, lwd = 15)

points(ndata$lxi[Levr], ndata$wi[Levr], col = "blue", pch = 25, lwd = 8)

points(ndata$lxi[StdRes], ndata$wi[StdRes], col = "green", pch = 20, lwd = 5)

text(ndata$lxi[as.numeric(row.names(outs))], ndata$wi[as.numeric(row.names(outs))], labels = round(ndata$wi[as.numeric(row.names(outs))], 3), pos = 1)

identify(ndata$lxi, ndata$wi)

and got the below, but not sure how I identify the points: enter image description here

my questions are:

  1. How to identify to deal with the outliers! After correcting for the leverage point error (all points far from leverage point)

  2. Do I need to use “plm” or it’s ok since I identified my data as panel?

  3. How to do two way tests for Random?

  4. Testing for heteroskedasticity by:

p<-plm(wi~lxi + Class+Class:lxi, data=ndata1)



the null was rejected(Null: homoskedastic). Does this solve the heteroskedasticity problem:

vcovHC(p, omega = Null, type="HC4")

coeftest(p, df=Inf, vcov=vcovHC(p, type="HC4"))

or I need to remove the insignificant points!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.