#Soal1
#PACKAGE
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(ggplot2)
library(car)
## Loading required package: carData
#MEMANGGIL DATA
data=read.table(file.choose(),header=T)
data
##    Petak    Y  X1  X2
## 1      1 34.4 100 140
## 2      2 38.7 110 180
## 3      3 34.2 120 150
## 4      4 41.7 130 200
## 5      5 39.0 140 170
## 6      6 41.5 150 220
## 7      7 40.9 160 190
## 8      8 44.0 170 240
## 9      9 43.9 180 200
## 10    10 47.6 190 230
## 11    11 49.7 200 269
## 12    12 48.7 210 220
#SCATTER PLOT Y (Hasil Panen Padi) DAN X1 (Dosis Pupuk)
ggplot(data, aes(x=X1, y=Y)) +
  geom_point(color="hotpink", size=2) +
  geom_smooth(method="lm", se=FALSE,
              color="springgreen4", linewidth=1.2) +
  ggtitle("Scatter Plot Antara Dosis Pupuk dan Hasil Panen Padi") +
  ylab("Hasil Panen Padi (kuintal/ha)") +
  xlab("Dosis Pupuk (kg/ha)") +
  theme(plot.title=element_text(hjust=0.5))
## `geom_smooth()` using formula = 'y ~ x'

#SCATTER PLOT Y (Hasil Panen Padi) DAN X2 (Curah Hujan)
ggplot(data, aes(x=X2, y=Y)) +
  geom_point(color="hotpink", size=2) +
  geom_smooth(method="lm", se=FALSE,
              color="springgreen4", linewidth=1.2) +
  ggtitle("Scatter Plot Antara Curah Hujan dan Hasil Panen Padi") +
  ylab("Hasil Panen Padi (kuintal/ha)") +
  xlab("Curah Hujan (mm/bulan)") +
  theme(plot.title=element_text(hjust=0.5))
## `geom_smooth()` using formula = 'y ~ x'

#REGRESI LINIER BERGANDA
model=lm(Y~X1+X2,data=data)
summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92924 -0.97284 -0.00117  1.03601  1.70057 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.39287    2.21841   7.840  2.6e-05 ***
## X1           0.07918    0.01942   4.077  0.00277 ** 
## X2           0.06157    0.01862   3.307  0.00913 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.356 on 9 degrees of freedom
## Multiple R-squared:  0.9415, Adjusted R-squared:  0.9285 
## F-statistic: 72.45 on 2 and 9 DF,  p-value: 2.828e-06
#PENDETEKSIAN MULTIKOLINIERITAS
#VIF < 10
vif(model)
##       X1       X2 
## 2.930714 2.930714
#PENGUJIAN ASUMSI
#Uji Heteroskedastisitas
bptest(model,varformula=~fitted.values(model),studentize=FALSE)
## 
##  Breusch-Pagan test
## 
## data:  model
## BP = 0.32528, df = 1, p-value = 0.5685
dwtest(model,alternative='two.sided')
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.283, p-value = 0.6963
## alternative hypothesis: true autocorrelation is not 0
shapiro.test(resid(model))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.9462, p-value = 0.5822
#Soal2
#Mengubah data menjadi matriks
n<-nrow(data)
n
## [1] 12
k<-ncol(data)
k
## [1] 4
#Membuat matriks
x0<-as.matrix(rep(1,n))
x1<-as.matrix(data$X1)
x2<-as.matrix(data$X2)
X<-as.matrix(cbind(x0,x1,x2))
X
##       [,1] [,2] [,3]
##  [1,]    1  100  140
##  [2,]    1  110  180
##  [3,]    1  120  150
##  [4,]    1  130  200
##  [5,]    1  140  170
##  [6,]    1  150  220
##  [7,]    1  160  190
##  [8,]    1  170  240
##  [9,]    1  180  200
## [10,]    1  190  230
## [11,]    1  200  269
## [12,]    1  210  220
Y<-as.matrix(data$Y)
Y
##       [,1]
##  [1,] 34.4
##  [2,] 38.7
##  [3,] 34.2
##  [4,] 41.7
##  [5,] 39.0
##  [6,] 41.5
##  [7,] 40.9
##  [8,] 44.0
##  [9,] 43.9
## [10,] 47.6
## [11,] 49.7
## [12,] 48.7
#Penaksiran parameter dengan matriks
x.trans<-t(X)
a<-x.trans %*% X
a
##      [,1]   [,2]   [,3]
## [1,]   12   1860   2409
## [2,] 1860 302600 385500
## [3,] 2409 385500 499161
a.inv<-solve(a)
a.inv
##               [,1]          [,2]          [,3]
## [1,]  2.6746125391  0.0002525866 -0.0131030144
## [2,]  0.0002525866  0.0002049450 -0.0001594972
## [3,] -0.0131030144 -0.0001594972  0.0001884188
b<-x.trans %*% Y
b
##          [,1]
## [1,]    504.3
## [2,]  80044.0
## [3,] 103154.3
b.topi<-a.inv %*% b
b.topi
##             [,1]
## [1,] 17.39287316
## [2,]  0.07917565
## [3,]  0.06156862
#Membuat nama koefisien
var_names<-c("intercept","Beta_1","Beta_2")

#Membuat tabel hasil
hasil<-data.frame(
  Estimate=as.numeric(b.topi),
  row.names = var_names
)
print(hasil)
##              Estimate
## intercept 17.39287316
## Beta_1     0.07917565
## Beta_2     0.06156862
#Membuat model regresi
model=lm(Y~X1+X2,data=data)
summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92924 -0.97284 -0.00117  1.03601  1.70057 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.39287    2.21841   7.840  2.6e-05 ***
## X1           0.07918    0.01942   4.077  0.00277 ** 
## X2           0.06157    0.01862   3.307  0.00913 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.356 on 9 degrees of freedom
## Multiple R-squared:  0.9415, Adjusted R-squared:  0.9285 
## F-statistic: 72.45 on 2 and 9 DF,  p-value: 2.828e-06
#MEMBENTUK MATRIKS
Y <- as.matrix(data$Y)
X <- as.matrix(cbind(data$X1,data$X2))
X <- cbind(1,X)
n <- nrow(X)
k <- ncol(X)-1

#KOEFISIEN DETERMINASI
Y.topi <- X %*% b.topi
Y.topi
##           [,1]
##  [1,] 33.93005
##  [2,] 37.18455
##  [3,] 36.12924
##  [4,] 39.99943
##  [5,] 38.94413
##  [6,] 42.81432
##  [7,] 41.75902
##  [8,] 45.62920
##  [9,] 43.95821
## [10,] 46.59703
## [11,] 49.78996
## [12,] 47.56486
JKT <- sum((Y-mean(Y))^2)
JKT
## [1] 283.1825
JKR <- sum((Y.topi-mean(Y))^2)
JKR
## [1] 266.6224
R2 <- JKR/JKT
R2
## [1] 0.9415214
#UJI F DENGAN MATRIKS
F_hitung <- (R2/k)/((1-R2)/(n-k-1)) 
F_hitung
## [1] 72.45123
#NILAI F TABEL
df1 <- k
df2 <- n-k-1
f_tabel <- qf(0.05,df1,df2,lower.tail = FALSE)
f_tabel
## [1] 4.256495