Først skal data indlæses og konveteres til R-data

# load pakken der kan indlæse statafiler (pakken hedder foreign)

library(foreign)
data<- read.dta(file.path("C:/Users/Vision/Desktop/STATISTIK/BIOSTATISTIK2/exercise20-1.dta"))
data
      Y   X1   X2   X3   X4
1  51.4  0.2 17.8 24.6 18.9
2  72.0  1.9 29.4 20.7  8.0
3  53.2  0.2 17.0 18.5 22.6
4  83.2 10.7 30.2 10.6  7.1
5  57.4  6.8 15.3  8.9 27.3
6  66.5 10.6 17.6 11.1 20.8
7  98.3  9.6 35.6 10.6  5.6
8  74.8  6.3 28.2  8.8 13.1
9  92.2 10.8 34.7 11.9  5.9
10 97.9  9.6 35.8 10.8  5.5
11 88.1 10.5 29.6 11.7  7.8
12 94.8 20.5 26.3  6.7 10.0
13 62.8  0.4 22.3 26.5 14.3
14 81.6  2.3 37.9 20.0  0.5

(a)

# Så laves en model hvor alle faktorer indgår som variable

#Modellen er Y=a + B1X1+B2X2+B3X3+B4X4

model.full <- lm(Y ~ X1 + X2 + X3 + X4, data=data)
summary(model.full) # show results

Call:
lm(formula = Y ~ X1 + X2 + X3 + X4, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1388 -1.3761 -0.7588  2.1733  4.0201 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -30.1369    37.5282  -0.803  0.44264   
X1            2.0699     0.4562   4.537  0.00141 **
X2            2.5816     0.7397   3.490  0.00683 **
X3            0.6360     0.4603   1.382  0.20039   
X4            1.1060     0.7648   1.446  0.18208   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.107 on 9 degrees of freedom
Multiple R-squared:  0.9757,    Adjusted R-squared:  0.9648 
F-statistic: 90.18 on 4 and 9 DF,  p-value: 2.954e-07

Så modellen kommer altså til at se således ud:

Y=-30.137+2.070X1+2.582X2+0.636X3+1.106X4

(b)

#H-Null=no significant multiple regression relationship

Fra ovenstående ses F=90.18 og P-value=2.95e-07 => reject H-Null

(c)

# Kan ses i første output:

For 1+2 (P hhv 0.0141 og 0.0068) kan H-Nul forkastes men ikke for 3+4 (P hhv 0.2004 og 0.1820)

(d)

# svarende har vi igen fra første output

Residual standard error = 3.107 on 9 degrees of freedom

Adjusted R-squared (coefficient of determination) = 0.9648

(e)

Y=-30.137 + 5.2*2.070 + 21.3*2.582 + 19.7*0.636 + 12.2*1.106

Y=61.646 (g)

Residual plot

# plot af residualer
plot(residuals(model.full))

#residualer vs fitted values
layout(matrix(c(1,2,3,4),2,2))
plot(model.full)

### Det ses at betingelserne for normalfordeling er opfyldt (øverste venstre plot samt nederste venstre plot) ### Det er svært helt at forkaste heteroscedasitet, jeg synes grafisk ikke det er til stede og derfor er betingelserne opfyldt. Der findes vist test der kan spytte et tal ud for dette.

Variable selection

#Backward variable selection

#Det ser ud til at X3 er mindst signifikant, så vi starter med at smide den ud af modellen
model.reduced1 <- lm(Y ~ X1 + X2 + X4, data=data)
summary(model.reduced1) # show results

Call:
lm(formula = Y ~ X1 + X2 + X4, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4107 -1.9354 -0.4523  2.5472  4.4035 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  16.8304    16.6111   1.013  0.33485    
X1            1.4769     0.1614   9.149 3.57e-06 ***
X2            1.7329     0.4306   4.024  0.00242 ** 
X4            0.2138     0.4281   0.499  0.62838    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.245 on 10 degrees of freedom
Multiple R-squared:  0.9705,    Adjusted R-squared:  0.9616 
F-statistic: 109.6 on 3 and 10 DF,  p-value: 5.981e-08
#Når man kigger på adjusted R-squared er det ikke den store forskel på de to modeller. Siden den antageligt ikke bidrager nævneværdigt med information til modellen smider vi den ud.
#Næste variabel der forsøges smidt ud af modellen er X4:
model.reduced2 <- lm(Y ~ X1 + X2, data=data)
summary(model.reduced2) # show results

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

Residuals:
   Min     1Q Median     3Q    Max 
-4.609 -1.510 -0.799  2.341  4.827 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  24.9632     3.1426   7.944 6.99e-06 ***
X1            1.4761     0.1558   9.473 1.27e-06 ***
X2            1.5264     0.1157  13.195 4.36e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.132 on 11 degrees of freedom
Multiple R-squared:  0.9698,    Adjusted R-squared:  0.9643 
F-statistic: 176.4 on 2 and 11 DF,  p-value: 4.4e-09
#Igen når man kigger på adjusted R-squared er det ikke den store forskel på de to modeller. Siden den antageligt ikke bidrager nævneværdigt med information til modellen smider vi den ud. VI kan se at P-værdien bliver mindre som konsekvens af at der er flere DF i de reducerede modeller.
#Så prøver vi at smide X2 ud - her må man forvente at det påvirker modellens egnethed:
model.reduced3 <- lm(Y ~ X1, data=data)
summary(model.reduced3) # show results

Call:
lm(formula = Y ~ X1, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-18.5817  -9.3153  -0.2446   7.5992  16.6881 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  62.3085     5.3632  11.618 6.93e-08 ***
X1            2.0108     0.5909   3.403  0.00524 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.3 on 12 degrees of freedom
Multiple R-squared:  0.4911,    Adjusted R-squared:  0.4487 
F-statistic: 11.58 on 1 and 12 DF,  p-value: 0.005241
# Ganske rigtigt falder Adjusted R-squared betragteligt (og P-værdien stiger markant) og X2 må holdes i modellen.

Vi ender med denne model Y=24.96+1.48X1+1.53X2 med 11 DF og P=4.4e-09