Linjär regression

En regressionsanalys används om man är intresserad av ett orsakssamband mellan två kontinuerliga variabler. Den man är intresserad av att förklara och som är y-axeln i en graf kallas responsvariabel (eller beroende variabel) medan den variabel som man vill förklara med kallas prediktor (eller oberoende variabel) och utgör x-axeln.

Den enklaste typen av regression är en linjär regression. När man behärskar den linjära regressionen kan den enkelt byggas ut för att även studera icke-linjära samband, olika regressionslinjer för olika grupper och ännu mer avancerade samband.

En linjär regression har följande krav på datat

Residual En residual är skillnaden mellan ett mätvärda av responsvariabeln och det förutsagda värdet av responsvariabeln om det hade varit placerat exakt på regressionslinjen. Ju längre från linjen (i y-led) en punkt ligger, desto större residual.

Läs in data

Vi vill undersöka reproduktionen hos en växt. Därför har vi mätt höjden på 20 plantor och vi har även räknat hur många frön varje planta har. Nu vill vi ta reda på om antalet frön beror på plantornas höjd.

Ladda ner följande fil plant_reproduction.txt (högerklicka, välj “spara länk som”) och spara filen på din hårddisk i en mapp med ett lämpligt namn.

Fortsätt med att läsa in datasetet och ge det ett namn, i det här fallet kallar vi det plant_reproduction_data. En detalerad beskrivning i hur man läser in filer finns i vår tidigare tutorial Läsa in data i R.

Glöm inte att dokumentera din kod i ett script, med kommentarer som förklarar vad du gör! Se vår tutorial om script om du behöver påminnelse om hur man skapar och använder script.

plant_reproduction_data <- read.table("data/plant_reproduction.txt",
                           header=T,
                           sep="\t",
                           dec=",")

<- betyder att vi sparar resultatet av funktionen (i det här fallet read.table() i ett objekt. Pilen pekar på objektet, och vi ger objektet ett namn, i det här fallet plant_reproduction_data.

read.table är kommandot för att läsa in textfiler med ändelsen .txt

"data/plant_reproduction.txt" är sökvägen till platsen där datafilen är sparad på din dator(data/), samt namnet på filen (plant_reproduction.txt).

header=Tbetyder att första raden i filen är kolumnernas rubriker, dvs inte värden.

sep="\t" betyder att alla värden är separerade med TAB.

dec="," betyder att komma används som decimalavgränsare (ex. 1,32).

Inspektera data

Börja med att titta på datans struktur med str().

Code 
str(plant_reproduction_data)
'data.frame':   20 obs. of  2 variables:
 $ Seeds       : int  466 501 467 489 496 521 522 517 506 427 ...
 $ Plant_height: int  45 45 44 54 52 57 53 50 48 32 ...

$ Seeds: int betyder att värderna i kolumnen Seeds är heltal (integers)

$ Plant_height: int betyder att värderna i kolumnen Mass är heltal (integers)

Om du har en kolumn som innehåller decimaltal skall det stå num. Om det istället står chr betyder det att datat har punkt som decimalavgränsare, medan du läste in det med komma som avgränsare (eller tvärtom) och R tror att det är text. Justera koden och läs in datat på nytt (se vår tidigare tutorial Läsa in data i R för justerad kod).

Visa sedan de fem första raderna av ditt dataset med head() för att se att allt ser korrekt ut

Code 
head(plant_reproduction_data)
  Seeds Plant_height
1   466           45
2   501           45
3   467           44
4   489           54
5   496           52
6   521           57

Visualisera datat med en enkel graf

För att göra grafen måste vi bestämma oss vad som skall vara på x-axeln, och vad som skall vara på y-axeln. Vilket av följande påstående är mest logiskt?

  1. Växtens storlek påverkar hur många frön den har

  2. Antalet frön en växt har påverkar dess storlek

Alternativ A är det som är mest logiskt. Alltså är det antalet frön som beror av växtens storlek. Därför skall vi ha antalet frön (Seeds) på y-axeln, och växtens storlek (Plants_height) på x-axeln.

Vi gör en enkel graf med funktionen plot() och lägger sedan till en regressionslinje med funktionen abline(). Grafen är inte publikationsduglig men ger en snabb överblick av våra data. Vi går igenom hur man gör en publiceringsduglig graf i den avslutande delen av vår tutorial, under Publiceringsduglig figur.

plot(Seeds ~ Plant_height, data = plant_reproduction_data)
abline(lm(Seeds ~ Plant_height, data = plant_reproduction_data))

Hur tolkar du grafen? Verkar det som att plantans höjd påverkar antalet frön den har?

Statistisk modellering

Vi vill nu göra en regressionsanalys för att undersöka om vår responsvariabel (beroende variabel) Seeds beror av vår förklarande variabel (oberoende variabel) Plants_height.

Vi specificerar en linjär modell med hjälp av funktionen lm() och väljer att spara resultatet i ett objekt som vi kallar m.plant_reproduction. Jag föredrar att alla mina modeller (resultat av statistka test) har ett namn som börjar med m. för att jag skall veta vad som är dataset och vad som är modeller. Ge alltid dina modeller beskrivande namn.

m.plant_reproduction<-lm(Seeds~Plant_height,
                    data=plant_reproduction_data)

I vår modell har vi vår responsvariabel (det vi har på y-axeln) till vänster om tilde-tecknet ~ och vår förklarande variabel till höger. data=plant_reproduction_data betyder att vi använder oss av datasetet vi tidigare döpt till plant_reproduction_data.

Statistiska resultat

Vi börjar med en ANOVA-tabell, och använder funktionen anova()

anova(m.plant_reproduction)
Analysis of Variance Table

Response: Seeds
             Df Sum Sq Mean Sq F value    Pr(>F)    
Plant_height  1  18361 18360.8  68.681 1.477e-07 ***
Residuals    18   4812   267.3                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vi får en ANOVA-tabell och kan inspektera resultatet. Vi ser att p-värdet (här kallat Pr(>F)) är mindre än 0.05, dvs Plant_height har en signifikant effekt på Seeds. R indikerar även att det är signifikant genom att lägga till en eller flera stjärnor *. När du rapporterar p-värden i skrift anger du p-värdet med tre decimaler, eller “p < 0.001” vid väldigt låga p-värden.

Är det en positiv eller negativ effekt av Plant_height påSeeds? För att veta det behöver vi ta fram skattningen av alla parametrar i modellen, dvs var x-axeln skär y-axeln, och linjens lutning. För det användervi oss av funktionen summary().

summary(m.plant_reproduction)

Call:
lm(formula = Seeds ~ Plant_height, data = plant_reproduction_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.376 -13.475   3.509  12.360  24.505 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  317.2300    22.2835  14.236 3.08e-11 ***
Plant_height   3.7064     0.4472   8.287 1.48e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 16.35 on 18 degrees of freedom
Multiple R-squared:  0.7923,    Adjusted R-squared:  0.7808 
F-statistic: 68.68 on 1 and 18 DF,  p-value: 1.477e-07

Coefficients ger oss fakta om koefficienterna. Jämför gärna med figuren ni gjorde ovan. (Intercept) är värdet på y-axeln när x-axeln är noll, dvs där regressionslinjen skär y-axeln. Vi ser att den skär y-axeln vid 317.2300 och att p-värdet (Pr(>|t|)) är 3.08e-11 det vill säga mindre än 0.05, dvs interceptet är signifikant skiljt från noll. Plant_height är lutningen på regressionslinjen. Den har lutningskoefficienten 3.7064 och eftersom den inte börjar med ett minustecken vet vi att det är ett positivt samband mellan Plants_height ochSeeds (plustecken skrivs aldrig ut)

R-square beskriver hur mycket av variationen i datat som förklaras av din statistiska modell. Kan gå från 0 (ingenting förklaras) till 1 (allt förklaras). R ger dig två värden, nämlingen Multiple R-squared och Adjusted R-squared. Använd Adjusted R-squared eftersom den tar hänsyn till antal parametrar i modellen.

Presentera din statistiska analys

Antalet frön ökad med ökad höjd hos växten (linjär regression, F = 68.681, d.f. = 1 och 18, p < 0.001, R2 = 0.78).

Man kan uttrycka den meningen på många sätt, ett alternativ är:

Ju högre växten är, desto fler frön har den (linjär regression, F = 68.681, d.f. = 1 och 18, p < 0.001, R2 = 0.78).

Utvärdera den statistiska modellen

Var modellen lämplig att använda för ditt dataset?

Vi utvärderar modellen genom diagnostiska grafer genom att använda funktionen plot() på vår statistiska modell.

plot(m.plant_reproduction)

Vi får fyra grafer att utvärdera, de två första är viktigast. Du kan behöva trycka upprepade gånger på ENTER för att se alla graferna (i R kommer de en och en).

Residuals vs Fitted bör visa en hyfsat rak linje. Den visar hur mycket residualerna (skillnaden i y-led mellan varje värde av din response och regressionslinjen) avviker från regressionslinjen. Residualerna motsvaras av cirklar i grafen. Om du har ett mönster i avvikelserna så betyder det att en rak regressionslinje inte var en lämplig statistisk modell.

Normal Q-Q visar om residualerna är normalfördelade. De bör följa den diagonala streckade linjen. Om de avviker på ett systematiskt sätt är residualerna inte perfekt normalfördelade, och vi kan behöva förändra modellen, exempelvis genom att transformera data.

Scale-Location illustrerar om variationen i datat är lika över alla värden. Om variationen ökar mycket åt höger (ett vanligt fall) så har vi större variation vid högre värden. Kan lösas genom att transformera data.

Residuals vs Leverage används för att hitta extremvärden som har onormalt stor påverkan på regressionslinjen. Mönstret i grafen är inte intressant, vi letar efter värden som ligger utanför de grå linjerna, speciellt linjerna för 1. Man bör dubbelkolla sådana värden (outliers) och fundera på om de skall vara med i datasetet. Kanske analysera såväl med som utan extremvärderna?

Publiceringsduglig figur

Vi avslutar med att göra en publiceringsduglig figur, och använder oss av paketet ggplot2. Om du inte sedan tidigare har paketet installerat så gör du det med koden install.packages("ggplot2"). Innan du använder paketet behöver du läsa in det i din session i R genom funktionen library()

Först - en enkel plot med ggplot2

library(ggplot2)

plot.plant_reproduction.simple <- ggplot(plant_reproduction_data, aes(x = Plant_height, y = Seeds)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic()

plot.plant_reproduction.simple
`geom_smooth()` using formula = 'y ~ x'

Vi går nu igenom koden för vår plot. Den är utspridd över flera rader, där varje rad avslutas med ett + för att R skall förstå att koden fortsätter på nästa rad.

ggplot(plant_reproduction_data, aes(x = Plant_height, y = Seeds)) + betyder att vi använder oss av funktionen ggplot(), och det första vi anger är namnet på vårt dataset, som vi kallade plant_reproduction_data.

aes() är en funktion inom ggplot, här anger vi estetiken (aes står för aestetics). Vi anger namnet på de variabler som skall vara på x-axeln respektive y-axeln med x = Plant_height, y = Seeds.

geom_point() + betuder att vi vill att varje datapunkt skall bli en punkt

geom_smooth(method = "lm", se = FALSE) + betyder att vi lägger till en linjär regressionslinje med method = "lm" och vi vill inte ha standardavvikelser, vilket vi anger med se = FALSE. Man brukar inte ha med det i grafer med regressioner. Vill du ändå ha osäkerheten i grafen kan du ersätta FALSE med TRUE.

theme_classic() är ett tema för grafer. Det är det temat som gör att graferna ser ut så som vi förväntar oss i naturvetenskapliga publikationer (vit bakgrund, inga stödlinjer, ren design). Notera att det är sista raden i vårt script, så den skall inte avslutas med plustecken.

Publiceringsbar figur

Den blev helt ok, men vi kan förbättra den.

  1. Vi vill att det skall stå “Plant height (cm)” på x-axeln, fixas med xlab()

  2. Vi vill att det skall stå “Number of seeds” på y-axeln, fixas med ylab()

  3. Vi vill välja blå punkter, och vi vill att de skall vara lite större,fixas med cex och color inom geom_point()

  4. Vi vill att regressionslinjen skall vara svart och som punkter, fixas med color och lty (står för line type) i geom_smooth()

plot.plant_reproduction <- ggplot(plant_reproduction_data, aes(x = Plant_height, y = Seeds)) +
  geom_point(cex = 3, colour = "blue") +
  geom_smooth(method = "lm", se = FALSE, lty = 3, color = "black") +
  xlab("Plant height (cm)")+
  ylab("Number of seeds")+
  theme_classic()

plot.plant_reproduction
`geom_smooth()` using formula = 'y ~ x'