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!

| cite | improve this question | | | | |

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.