Station Charge

list.files()

library(readxl)


station_data_fix <- read_excel("C:/Users/Admin/Downloads/station_data_fix.xlsx")


station_data_fix = as.data.frame(station_data_fix)


station_data_fix

NA

Y<-station_data_fix$`dollars`
X1<-station_data_fix$`kwhTotal`
X2<-station_data_fix$`chargeTimeHrs`

data1 <- data.frame(Y,X1,X2)

#Uji Multikolinearitas
library(car)
vif(lm(Y~X1+X2, data=data1))

      X1       X2 
1.018421 1.018421 

#Model Poisson
model_poisson <- glm(formula=Y~X1+X2, family=poisson(link="log"), data=data1)
summary(model_poisson)

Call:
glm(formula = Y ~ X1 + X2, family = poisson(link = "log"), data = data1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2262  -0.6042  -0.5826  -0.5304   2.3845  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.07081    0.30473  -6.796 1.08e-11 ***
X1           0.03136    0.04013   0.781    0.435    
X2           0.05361    0.01183   4.530 5.90e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 179.28  on 249  degrees of freedom
Residual deviance: 168.30  on 247  degrees of freedom
AIC: 257.21

Number of Fisher Scoring iterations: 6

library(AER)
dispersiontest(model_poisson)

#ZIP
library(pscl)
model_zip <- zeroinfl(formula=Y~X1+X2, data=data1)

Warning: glm.fit: algorithm did not convergeWarning: glm.fit: fitted probabilities numerically 0 or 1 occurred

summary(model_zip)

Call:
zeroinfl(formula = Y ~ X1 + X2, data = data1)

Pearson residuals:
       Min         1Q     Median         3Q        Max 
-3.573e-01 -1.120e-08 -1.100e-08 -1.083e-08  7.193e-01 

Count model coefficients (poisson with log link):
              Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.330e-02  5.205e-01  -0.026    0.980
X1           3.234e-02  7.526e-02   0.430    0.667
X2           1.108e-05  1.904e-02   0.001    1.000

Zero-inflation model coefficients (binomial with logit link):
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   6742.861 342760.102   0.020    0.984
X1              -3.787   1632.564  -0.002    0.998
X2           -1682.611  84061.990  -0.020    0.984

Number of iterations in BFGS optimization: 11 
Log-likelihood: -43.68 on 6 Df

library(VGAM)

#Membangkitkan Respon
p<-predict(model_zip,type="zero")
lambda<-predict(model_zip,type="count")
koef<-model_zip$coefficients[2]
beta0<-koef$zero[1]
beta1<-koef$zero[2]
beta2<-koef$zero[3]

#Simulasikan
simul <-function(X1,X2,Y) {
  model_zip1<-zeroinfl(Y~X1+X2)
  beta<-model_zip1$coefficients[2]
  beta_z<-beta$zero
}

iterasi <- 600
betaduga0<-rep(0,iterasi)
betaduga1<-rep(0,iterasi)
betaduga2<-rep(0,iterasi)

for (i in 1:iterasi) {
  Y<-rzipois(250,lambda = lambda, pstr0 = p)
  betaduga0[i] <-simul(X1,X2,Y)[1]
  betaduga1[i] <-simul(X1,X2,Y)[2]
  betaduga2[i] <-simul(X1,X2,Y)[3]
}

#Mentaksir Koefisien
betaduga<-data.frame(betaduga0, betaduga1, betaduga2)
names(betaduga)<-c("beta0","beta1","beta2")
head(betaduga)

#Plot (Hasil Simulasi)
plot(betaduga0, main=expression(paste("Hasil penaksiran ", beta[0])),
     ylab=expression(beta[0]), xlab="Iterasi")

plot(betaduga1, main=expression(paste("Hasil penaksiran ", beta[1])),
     ylab=expression(beta[1]), xlab="Iterasi")

plot(betaduga2, main=expression(paste("Hasil penaksiran ", beta[2])),
     ylab=expression(beta[2]), xlab="Iterasi")