Statystyka w R

Oprócz analizy rozkładu jednej zmiennej często w badaniach statystycznych rozważany jest problem analizy zależności pomiędzy zmiennymi. Analiza zależności obejmuje m.in. analizę korelacji i wyznaczanie funkcji regresji.

####################################################################
############## Analiza zależności - wybrane zagadnienia ############
####################################################################

### wspolczynnik korelacji liniowej Pearsona - wersja 1 programu
wsp_kor=function(x,y) {
n=length(x)
xsr=mean(x)
ysr=mean(y)
kowariancja=sum((x-xsr)*(y-ysr))/n
wariancja_x=sum((x-xsr)^2)/n
wariancja_y=sum((y-ysr)^2)/n
odch_stand_x=sqrt(wariancja_x)
odch_stand_y=sqrt(wariancja_y)
r=kowariancja/(odch_stand_x*odch_stand_y)
r
}

### wspolczynnik korelacji liniowej Pearsona - wersja 2 programu
wsp_kor=function(x,y) {
n=length(x)
kowariancja=cov(x,y)*(n-1)/n
wariancja_x=var(x)*(n-1)/n
wariancja_y=var(y)*(n-1)/n
odch_stand_x=sqrt(wariancja_x)
odch_stand_y=sqrt(wariancja_y)
r=kowariancja/(odch_stand_x*odch_stand_y)
r
}

#wykorzystujemy zbiór women
attach(women) # nazwy kolumn jako nazwy zmiennych
wsp_kor(weight,height)
# otrzymujemy:
# [1] 0.9954948
# lub korzystając z gotowej funkcji w R, tj.:
cor(weight,height)

### Szacowanie parametrów funkcji regresji
### z jedną zmienną objaśniającą i ze stałą
funkcja_reg=function(x,y){
n=length(x)
xsr=mean(x)
ysr=mean(y)
kowariancja=sum((x-xsr)*(y-ysr))/n
wariancja_x=sum((x-xsr)^2)/n
a1=kowariancja/wariancja_x # ocena parametru kierunkowego
a0=ysr-a1*xsr # ocena stałej
yteor=a1*x+a0 # ciag wartości teoretycznych
Su2=sum((y-yteor)^2)/(n-2) # wariancja resztowa
Su=sqrt(Su2) #odchylenie standardowe reszt
R2=1-(sum((y-yteor)^2)/sum((y-ysr)^2)) # współczynnik determinacji
#wyniki jako lista
list(a1=a1,a0=a0,Su=Su,R2=R2,yteor=yteor)
}
# Przykład:
attach(women)
funkcja_reg(weight,height)
# otrzymujemy:
# $a1
# [1] 0.2872492
# $a0
# [1] 25.72346
# $Su
# [1] 0.4400391
# $R2
# [1] 0.9910098
# $yteor
# [1] 58.75712 59.33162 60.19336 61.05511 61.91686 62.77861 63.64035 64.50210
# [9] 65.65110 66.51285 67.66184 68.81084 69.95984 71.39608 72.83233

#te same rezultaty można uzyskać następująco:
lm(height~weight)
#lub dokładniej:
summary(lm(height~weight))
#lub wybrane wyniki:
summary(lm(height~weight))$coef # współczynniki (estymatory parametrów)
summary(lm(height~weight))$coef[,1] # pierwsza kolumna
summary(lm(height~weight))$sigma # odchylenie standardowe reszt
summary(lm(height~weight))$r.squared # współczynnik determinacji
#wykaz dostępnych elementów listy wyników
attributes(summary(lm(height~weight)))
#ponadto:
fitted(lm(height~weight)) # wartości teoretyczne