Multivariata metoder - Uppsala Universitet
Jacob Lindberg
20 Oct 2015
Detta dokument finns tillgängligt på http://rpubs.com/lindberg/118713
- INLEDNING
- IMPORTERA DATA
- FRÅGA A
- FRÅGA B - Description of any changes made to the dataset, - Short description of the statistical technique(s) - Formulation of model - Formulation of hypotheses - Presentation of results - Conclusions
- FRÅGA C
- FRÅGA D
INLEDNING
IMPORTERA DATA
Steg 1 är att exportera från SAS till CSV.
Steg 2 är att ladda in våra dataset.
df <- read.csv("C:/Users/jl/Dropbox/aktuellaKurser/multivmetoder/exammm/homedata.csv")Dessa paket kommer att användas.
library(car) #for regression diagnostics
library(dplyr) #for data minupulation##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionFRÅGA A
MÅL
The municipality is interested in the following:
- The relationship between knowledge variables and kindergarten experience, taking background variables into account.
- Is there any grouping of individuals that have certain combinations of variables in common, and if so: what is it that distinguishes these groups?
Vårt mål är att se hur resultat på math, literacy, social samvarierar med intressevariabler samtidigt som vi håller bakgrundsvariabler konstant.
BESKRIVNING AV VARIABLER
Våra variabelnamn är
names(df)## [1] "id" "birth_year" "birth_month"
## [4] "kindergarten_type" "kindergarten_amt" "literacy"
## [7] "social" "math" "parents_educ"
## [10] "gender"Vi delar in dessa i
- Y-variabler (dependent) är
math, literacy, social. De är alla resultat på test. - X-variabler (explanatory) delas in i
- Intressevariabler:
kindergarten_type, kindergarten_amt. - Bakgrundsvariabler:
gender, parents_educ, birth_year, birth_month, ID.
- Intressevariabler:
Anledningen för indelningen följer av frågans formulering “the relationship between the knowlegde variables (Y) and kindergarten experience (intressevars), taking background variables (bakgrundvars) into account.”
Nedan ser vi tabulering över kindergarten_type som är 1 om eleven går i Montessori, 2 om eleven går i Ur och skur, 3 om eleven går vanlig skola eller inte går i kindergarten.
table(df$kindergarten_type)##
## 1 2 3
## 38 6 452(Notera att vi ombads koda om till 1 2 3 4 i uppgiften, men då vi senare använder en dummy-kodning utesluts detta steg - det vore att göra en sak som ej används. Och vi demontrerar redan kunskapen att koda om variabler så det behövs ej.)
Nedan ser vi tabulering över kindergarten_amt som är 1 om eleven inte går på dagis, 2 om den går 15h, 3 om den går mer än 15h. De flesta elever går mer än 15 timmar. 1
table(df$kindergarten_amt)##
## 1 2 3
## 48 179 269Nedan ser vi en tabell över när eleverna är födda. (Se kod-kommentar för vilket fel som uppstår. Kod-kommentarer är att betrakta som en del i inlämningen.)
table(df$birth_year) #ser ett error i 20009. #ser att 2 pers har gått om. de är alltså speciella observationer och skulle därmed kunna vara outliers att ta bort, ifall vi vill göra inferens om elever som ej gått om. ##
## 2009 2010 20009
## 493 2 1table(df$birth_month) ##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 60 52 43 52 39 95 46 35 19 21 20 14Nedan er vi en tabell över gender där 2 är en pojke och 1 en flicka. Vi ser också en tabell över parents_educ som är antalet föräldrar med universitetsubildning.
table(df$gender)##
## 1 2
## 244 252table(df$parents_educ)##
## 0 1 2
## 227 214 55SIGNIFIKANS-NIVÅ
Vi väljer signifikans-nivå 5 procent i alla våra test, för att det är standard i social science.
QUICK EXPLORATORY DATA ANALYSIS
Vi börjar med att - för våra Y - kolla range, rita boxplot, histogram, quantile quantile plots, samt kolla om missing values finns.
c(
range(df$math) ,
range(df$literacy) ,
range(df$social)
)## [1] 12 53 2 45 14 66#boxplot
#Dessa innhåller en-två outliers i varje (kan vara samma person som upprepas). Men de ligger mycket nära övriga elever. De ska ej uteslutas då den nog är korrekt rapporterad, det är blott en-två duktig(a) elev i ämnet i fråga.
par(mfrow=c(1,3))
boxplot(df$literacy, main="BOXPLOTS Literacy")
boxplot(df$social, main="Social")
boxplot(df$math, main="Math")
#histogram (mer detalj än boxplot) samt q q plot.
par(mfrow=c(3,2))
hist(df$literacy, main="BOXPLOTS - Literacy")
qqPlot(df$literacy, main="QUANTILE QUANTILE PLOTS - Literacy")
hist(df$social, main="Social")
qqPlot(df$social, main="Social")
hist(df$math, main="Math")
qqPlot(df$social, main="Math")
# alla tre ser NF ut.
par(mfrow=c(1,1)) #återställer mfrow
sum(is.na(df)) #==0 inga missing values## [1] 0EDITERA VARIABLERNA
Vi ändrar lite i variablerna så att vi
- kastar ut 2010 observationerna. Vi anser att något är speciellt med varför de är yngre. Vi kan slänga in dem i kvartal4, men detta känns mer korekt.
- skapar boy=1 om pojke, =0 om flicka.
- tar bort felkodning år 20009.
# year 2010 throw out ####
nrow(subset(df, df$birth_year == 2010)) #antal rader i det subset som har birth year == 2010## [1] 2df <- df[df$birth_year==2009,] #3 obs droppas
# BOY ####
#skapa boy=1 if gender=2
#syfte: skönare tolkning i regressionen
df$boy <- NULL #skapa
df$boy[df$gender == 2] <- 1
df$boy[df$gender == 1] = 0
#så boy=1 if pojke =0 if flicka.
table(df$boy, df$gender) #checking##
## 1 2
## 0 241 0
## 1 0 252# 2009 ####
#tar bort 20009 flekodning
df$birth_year[df$birth_year == 20009] <- 2009
table(df$birth_year) #checking##
## 2009
## 493ARGUMENTATION FÖR KVARTAL
Teorin säger att de som är födda tidigt på året presterar bättre. Vi väljer att dra gräns för “tidigt” som kvartal. Vi låter baseline vara fjärde kvartalet och kan på så sätt se om tex kvartal 1 har större coef än kvartal 2 (som är sign.) ty då kan vi säga att kvartal 1 presterar bättre.
#skapa vars
df$kvartal1 <- NULL
df$kvartal2 <- NULL
df$kvartal3 <- NULL
df$kva1 <- ifelse( df$birth_month < 4, 1, 0) #får 1 2 3
df$kva2 <- ifelse(df$birth_month > 3 & df$birth_month < 7, 1, 0) #får 4 5 6
df$kva3 <- ifelse(df$birth_month > 6 & df$birth_month < 10, 1, 0) #får 7 8 9
#kva4 = baseline
names(df)## [1] "id" "birth_year" "birth_month"
## [4] "kindergarten_type" "kindergarten_amt" "literacy"
## [7] "social" "math" "parents_educ"
## [10] "gender" "boy" "kva1"
## [13] "kva2" "kva3"#för pedagogikens skull visar vi nu hur kva ser ut som dummy.
df[1:10, c(3, 12:14)] #checking. denna väljer ut månad och alla kvartalen för 10 rader ## birth_month kva1 kva2 kva3
## 1 6 0 1 0
## 2 5 0 1 0
## 3 5 0 1 0
## 4 6 0 1 0
## 5 6 0 1 0
## 6 12 0 0 0
## 7 10 0 0 0
## 8 6 0 1 0
## 9 12 0 0 0
## 10 11 0 0 0ARGUMENTATION FÖR DUMMY
Notera names(df) output ovan, samt variable description. Vi ska nu fundera på vad som egentligen efterfrågas.
Vi vill mäta:
- om du går i kindergarten
- hur mycket du går i kinderkarten
- vilken sorts kindergarten
Det är våra intresservariabler. Hur påverkar de y1 y2 y3? Studera nu vars description på _amt och _type. Inger bad oss skapa =4 if didnt go to kindergarten, och det gjorde vi. Men allt det är onödigt… ty vi inför multikoll i modellen om vi gör så. Vi har då två ordinala variabler som försöker mäta tre saker, det är inte bra definierade variabler.
Notera alltså att vi går från två kateogivariabler till dummies, men bibehåller frågeställningen.
Vi skapar nu variablerna
- dagis = 1 om du gick på dagis, =0 om du inte gick på dagis.
- specialdagis = 1 om du gick på specialdagiset Montessori eller Iurochskur, =0 om du inte gick på specialdagis (dvs antingen vanligt dagis eller inget dagis).
- tidextra = 1 om du gick mer än 15h på dagis, =0 annars.
#dagis
df$dagis <- NULL
#dagis=1 om du går på dagis
df$dagis <- ifelse(df$kindergarten_amt == 1, 0, 1) #_amt=1 är didnt go to kindergarten.
table(df$dagis, df$kindergarten_amt)##
## 1 2 3
## 0 47 0 0
## 1 0 178 268#endast 47 pers som ej går på dagis
#TIDEXTRA? GIVET ATT DU GICK, GICK DU EXTRA LÄNGE?
df$tidextra <- NULL
#tidextra = 0 if 15h och =1 if över 15h.
df$tidextra <- ifelse(df$kindergarten_amt == 3, 1, 0)
#(notera att nu har både icke kinderfarten kids och de som gick i kindergarten med "bara" 15h =0 så denna variabel är en dummy som säger huruvida du gick extra länge när du var på kindergarten.)
table(df$tidextra, df$kindergarten_amt) #checking##
## 1 2 3
## 0 47 178 0
## 1 0 0 268#SPECIALDAGIS
#GIVET ATT DU GICK, GICK DU PÅ ETT SPECIALDAGIS?
df$specialdagis <- NULL
df$specialdagis <- ifelse(df$kindergarten_type == 3, 0, 1)
table(df$specialdagis, df$kindergarten_type) #checking. nr3 är other. nr1 är monte nr2 är uroskur. ##
## 1 2 3
## 0 0 0 449
## 1 38 6 0table(df$dagis, df$specialdagis) #3st är felaktigt inlagda i datan. 3 pers hade både dagis==0 och specialdagis==1.##
## 0 1
## 0 44 3
## 1 405 41df[df$dagis==0 & df$specialdagis==1, ] #visar sig att alla dessa kommer från montessori då _type=1## id birth_year birth_month kindergarten_type kindergarten_amt literacy
## 20 20 2009 4 1 1 18
## 33 33 2009 7 1 1 15
## 40 40 2009 10 1 1 18
## social math parents_educ gender boy kva1 kva2 kva3 dagis tidextra
## 20 30 13 0 2 1 0 1 0 0 0
## 33 19 13 0 2 1 0 0 1 0 0
## 40 14 20 0 2 1 0 0 0 0 0
## specialdagis
## 20 1
## 33 1
## 40 1#vi löser detta genom att sätta dem i specialdagis-gruppen. man kan antingen göra det eller sätta dem i vanliga gruppen.
df[df$dagis==0 & df$specialdagis==1, ]$dagis <- 1
table(df$dagis, df$specialdagis)##
## 0 1
## 0 44 0
## 1 405 44QUICK EXPLORATORY DEL 2
Först: våra Y.
Hur är våra Y korrelerade? Om de är det så är samvariationen stor, vilket kan vara intressant vid steget då vi ska tolka vår output.
pairs(~ math + literacy + social
#+ birth_month
, df)
#ca .20 i corr.
cor(df$literacy, df$social)## [1] 0.3072068cor(df$literacy, df$math)## [1] 0.2469654cor(df$math, df$social)## [1] 0.2097617Alltså är våra Y lite korrelerade med varandra.
Sedan: våra X.
table(df$tidextra)##
## 0 1
## 225 268table(df$dagis)##
## 0 1
## 44 449table(df$specialdagis)##
## 0 1
## 449 44table(df$birth_month) #denna ligger till grund för kvartal. för att se hur kvartal är uppbyggd se koden "df[1:10, c(3, 12:14)]"##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 58 51 43 52 39 95 46 35 19 21 20 14table(df$parents_educ)##
## 0 1 2
## 224 214 55sum(df$boy)## [1] 252Här ser vi att grupperna för våra förklarande variabler är för ojämna för att kunna göra Manova.
MOTIVERA VAL AV METOD
Uppgiften lyder “The relationship between knowledge variables and kindergarten experience, taking background variables into account.”
Datan vi har är
- Vi har Y1 Y2 Y3: literacy math social. Alla är numerical.
- intressevars = dagis, tidextra, specialdagis. Tre dummy.
- bakgrundvars = parents _ educ, boy. En kategorisk, en dummy.
Frågor att ställa oss innan vi bestämmer teknik (vi använder Ingers beslutsträd Enligt Lecture 1 slide 42):
- What type of relationships is being examined? Dependent or independent? Vi tycker dependent, för de tre Y vi har är sådant man typiskt sett vill se förklaras utav bakgrund. SVAR: DEPENDENT.
- Number of dependent. Vi har Y1 Y2 Y3 så SVAR: SEVERAL DEPENDENT IN SINGLE RELATIONSHIP.
- Measurement scale of the dependent variables. Våra Y går från 0 til 50. SVAR: NUMERICAL.
- Measurement scale of the explanatora variables. SVAR: Numerical -> cancorr och mutlivariate regression, SVAR: categorical -> manova. Vi har både numerical och categorical.
Då gruppstorlekarna är för små för att göra Manova (se kommentar tidigare också) samt att frågeställningen är “taking background variables into account” vilket regression hanterar mycket bra, så väljer vi att göra en multivariat multipel regression.
Det är viktigt att vi gör en multivariate multipel regression och inte bara en multipel regression, då “multivariate models are often more powerful, when the response variables are correlated.” Och vi ser av pairs() och cor() ovan (i QUICK EXPLORATORY DEL 2) att så är fallet för vår data.
SKAPA MODELLEN
Vi är intresserade av att tolka tid och specialdagis, medan vi håller badagisrunden kosntant
Nedan skapar vi vår multivariata multipla regression. \[math + literacy + social = dagis + specialdagis + tidextra + parentseduc + boy + kva1 +kva2 +kva3\]
Där cbind() gör en multivariate mulitple regression, i stället för en vanlig multiple regression, enligt http://www.psych.yorku.ca/lab/psy6140/lectures/MultivariateRegression2x2.pdf
fit1 <- lm(cbind(math, literacy, social) ~
dagis + specialdagis + tidextra #intressevariabler
+ parents_educ + boy + kva1 +kva2 +kva3 #kontrollvars
, data=df)
fit1.res <- residuals(fit1) #spara residuals(Vi gjorde hela analysen även med att inkludera interaktioner, men dessa var ej signifikanta.)
Finns ingen anledning att göra konfidensintervall då våra förklarande variabler är dummy.
Först kollar vi assumptions, sedan tolkar vi resultatet om assumptions är uppfyllda.
KOLLA ASSUMPTIONS
Sample size
Man ska gärna ha 20 gånger (antal förklaringsvariabler). Vid 10 förklaringsvariabler (vi har färre) så är n=200 rek. Vi har n=400 så sample size är bra.
Multikoll
Inger kollar tolerance typ II value. Vi kollar multikoll i stället - det samma sak egentligen då tolerance = 1 - r^2.
Korrelation med kategorisk ska ha Spearman. Vi har dummy så vi har Spearman, det är lite vanskligt att ha korrelation på dummy men dummy är vad vi har.
df_corrmatrix <- select(df, math, literacy, social, tidextra, dagis, specialdagis, boy, kva1, kva2, kva3)
cor(df_corrmatrix, method="spearman")## math literacy social tidextra dagis
## math 1.00000000 0.24328729 0.16831331 0.34232407 0.37734182
## literacy 0.24328729 1.00000000 0.29808839 0.55471348 0.28957799
## social 0.16831331 0.29808839 1.00000000 0.43256690 0.36129230
## tidextra 0.34232407 0.55471348 0.43256690 1.00000000 0.34164833
## dagis 0.37734182 0.28957799 0.36129230 0.34164833 1.00000000
## specialdagis 0.05549292 -0.04928029 -0.02445963 0.08686073 0.09799555
## boy -0.15897942 -0.17683134 -0.19493653 -0.10582393 -0.22073744
## kva1 0.46731458 0.33090144 0.27391463 0.36482733 0.19359441
## kva2 -0.19500624 -0.17277184 -0.01198926 -0.14376470 -0.06457797
## kva3 -0.20715609 -0.11695242 -0.21278723 -0.14543814 -0.10748579
## specialdagis boy kva1 kva2 kva3
## math 0.05549292 -0.15897942 0.46731458 -0.19500624 -0.20715609
## literacy -0.04928029 -0.17683134 0.33090144 -0.17277184 -0.11695242
## social -0.02445963 -0.19493653 0.27391463 -0.01198926 -0.21278723
## tidextra 0.08686073 -0.10582393 0.36482733 -0.14376470 -0.14543814
## dagis 0.09799555 -0.22073744 0.19359441 -0.06457797 -0.10748579
## specialdagis 1.00000000 -0.02121920 0.05290643 0.03522164 -0.10483005
## boy -0.02121920 1.00000000 -0.05883527 -0.03411501 0.08965387
## kva1 0.05290643 -0.05883527 1.00000000 -0.51967508 -0.33678158
## kva2 0.03522164 -0.03411501 -0.51967508 1.00000000 -0.39263681
## kva3 -0.10483005 0.08965387 -0.33678158 -0.39263681 1.00000000(1- 0.55471348^2) #tolerance lugnt på även den största corren. ## [1] 0.692293Slutsats: Ingen multikoll.
De fyra assumptions i linjär regressoin
- Linjäritet
- Homosked
- Oberoende (okorrelerade) feltermer
- Multivariat normalfördelade residualer
Dessa kollas nedan.
1 Linjäritet
“the true relationship is linear”
#skapa partial models (detta får sas automatiskt)
fit.math <- lm(math ~
dagis + specialdagis + tidextra
+ parents_educ + boy + kva1 +kva2 +kva3
, data=df)
fit.literacy <- lm(literacy ~
dagis + specialdagis + tidextra
+ parents_educ + boy + kva1 +kva2 +kva3
, data=df)
fit.social <- lm(social ~
dagis + specialdagis + tidextra
+ parents_educ + boy + kva1 +kva2 +kva3
, data=df)#skapa partial resid
res.math <- resid(fit.math)
res.literacy <- resid(fit.literacy)
res.social <- resid(fit.social)#plotta
#L10s34-36. vi letar efter non-linear patterns.
leveragePlots(fit.math, main="Leverage plots Math")
leveragePlots(fit.literacy, main="Leverage plots Literacy")
leveragePlots(fit.social, main="Leverage plots Social")
Av leveragePlots ser vi inga icke-linjära samband.
Slutsats: \(\vec{Y}\) är en linjär funktion av \(\vec{X}\) plus error term.
2 Homosked
“the error terms are random variables with mean 0 and constant variance (homosked)”
Vi testar detta grafiskt med hjälp av scatter på residualerna och quantile-quantile-plots.
#Samma plots som Inger L10s44
#modellnamn$fitted ger fitted values
#MATH PLOTS
par(mfrow=c(1,2))
plot(fit.math$fitted, fit.math$residuals, main="Math") #ser homosked ut
qqPlot(fit.math) #ser NF ut
#hist(fit.math$residuals) #ser NF ut
#LITERACY PLOTS
par(mfrow=c(1,2))
plot(fit.literacy$fitted, fit.literacy$residuals, main="Literacy") #ser homosked ut
qqPlot(fit.literacy) #ser NF ut
#hist(fit.literacy$residuals) #ser NF ut
#SOCIAL PLOTS
par(mfrow=c(1,2))
plot(fit.social$fitted, fit.social$residuals, main="Social") #ser homosked ut med tendens till
qqPlot(fit.social)
#hist(fit.social$residuals) #ser NF men tendens till lite skew
par(mfrow=c(1,1)) #återställMan skulle kunna se en liten tendens till heterosked i Literacy, speciellt i att stora värden på fitted ger negativa residuals. Man skulle kunna se en tendens till heterosked i Social, speciellt i att stora värden på fitted ger negativa residuals vilket även syns i quantile-quantile-plotten.
(Vi prövade att köra log(social) och log(math) och log(literacy) men det blev inte bättre.)
Slutsats: Inga allvarliga tecken på heterosked. Bra!
3 Oberoende feltermer
“the error term are note correlated, in other cords indep”
Detta har vi redan tittat i residualplot ovan. Vi ser inga allvarliga tecken på att feltermerna samvarierar.
4 Multivariat normalfördelade residualer
“the error terms are multivariate normally distributed”
Det är större slh att så är fallet om variablerna är NF (källa: L10 part 2 minut 14). Alla univariate är NF, så det är ingen överraskning att vår residualerna för multimodellen är NF enligt plot nedan.
qqPlot(fit1$resid, main="Quantile-quantile plot multimodell residuals")
Slutsats: Residualerna i fråga är normalfördelade.
FRÅGA B
Fråga A och B går in i varandra, och det är svårt att dela upp dem strikt. En del i A är relevant för B
Description of any changes made to the dataset,
e.g. deletion of any obvious errors, transformation of any variables, etc.
Vi började med en datacleaning, beskriven tidigare. Inklusive förändring av variabler och motivaton för detta. Utöver det har inga ändringar gjorts i denna del.
Vi fittade vår modell med log() på Y men det gav ej bättre uppfyllning av asusmptions.
Short description of the statistical technique(s)
I multivariat multipel linjär regression (MultiMLR) har vi regression på formen \[\vec{Y_i} = B \vex{X_i} + \vec{\epsilon_i}\] där \(B\) är en matris av våra beta.
I MultiMLR där vi kan ha fler än en Y-variabel. I vårt fall samvarierar våra Y, och det hade varit synd att bara göra tre MLR i stället för att göra tre stycket MultiMLR. Ty om man gör tre separa får följande negativa effekter:
- du kan inte ta hänsyn till korrelation mellan Y
- du kan inte göra simultana test - exempelvis kaske pojke=1 eller =0 inte spelar någon roll för varje univariat regression men i MultiMLR kan vi se att den faktiskt spelar roll.
- du kan inte göra hypotestestet \(H_0: \beta_{1i} = 0\) ensidigt, men det kan man i MultiMLR.
Formulation of model
(if appropriate, depending on chosen technique), or description of the investigated relationship(s) in a way that reflects the selected statistical technique
Vi sätter upp modellen \[Y_1 + Y_2 + Y_3 = intresse + bakgrund\] där \(bakgrund = parentseduc + boy + kva1 +kva2 +kva3\) och \(intresse = dagis + specialdagis + tidextra\) där dagis = 1 om du går till dags. specialdagis=1 om specialdagis. tidextra=1 om du gick 15+ timmar annars 0. Samt kva1, kva2, kva3 är en dummy för vilket kvartal där fjärde kvartalet är baseline. Vidare är boy=1 om pojke, och parentseduc mäter antalet föräldrar med universitetsutbildning.
Formulation of hypotheses
(if appropriate, depending on chosen technique)
Mycket viktigt att först titta på overall effect. Har vi totalt sett förklarande variabler? Sedan, och endast sedan, kan vi titta på resultatet på univariat nivå. Där t-test m.m. görs.
I Manova(fit1) nedan så görs F-testet \(H_0: \beta_i = 0\) mot \(H_1: någon \, \beta_i \ne 0\). Det p-värde som rapporteras är som vanligt, sannolikheten att observera ett lika extremt eller mer extremt värde givte att nollhypotesen är sann. Så vi vill att p-värdet ska vara litet för då finns det en effekt.
Library(car) paketets funktion Manova() använder Typ2 medan manova() i baseR använder Typ1. Källa: http://stats.stackexchange.com/questions/11127/multivariate-multiple-regression-in-r
När vi sedan konstaterat att en variabel är sign. på multinivå, då och endast då kan vi studera den på singelnivå (univariat nivå). I summary(fit1) nedan så görs t-testet \(H_0: t=0\) mot \(H_1: t \ne 0\). Det p-värde som rapporteras är som vanligt, sannolikheten att observera ett lika extremt eller mer extremt värde givte att nollhypotesen är sann. Så vi vill att p-värdet ska vara litet för då finns det en effekt.
Presentation of results
Har variablerna en effekt sign skiljd från noll? Ja.
#gör nu denna är ingers Lec10 slide 19, 21-22.
#notera att vi ej gör en manova-analys - uta nortfarande en Multi MReg - men använder Manova() för att skapa den output SAS får automatiskt.
Manova(fit1, test.statistic="Wilks")##
## Type II MANOVA Tests: Wilks test statistic
## Df test stat approx F num Df den Df Pr(>F)
## dagis 1 0.83680 31.335 3 482 < 2.2e-16 ***
## specialdagis 1 0.96800 5.311 3 482 0.001311 **
## tidextra 1 0.74566 54.802 3 482 < 2.2e-16 ***
## parents_educ 1 0.95158 8.176 3 482 2.561e-05 ***
## boy 1 0.96726 5.437 3 482 0.001103 **
## kva1 1 0.91398 15.120 3 482 2.027e-09 ***
## kva2 1 0.99193 1.307 3 482 0.271476
## kva3 1 0.99597 0.650 3 482 0.583559
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#längre output ges utav
#summary(Manova(fit1), test="Wilks") Av Manova(fit1) ser vi attde som är sign i multimodellen är dagis specialdagis tidextra parents _ educ boy kva1, där de som vi är intresserade av att tolka är dagis, specialdagis, tidextra dvs samtliga.
Detta är viktigt. Samtliga intressevars är sign. i multimodellen, nu kan vi fråga oss huruvida de är sign. i singelmodellerna (dvs på univariat nivå). Man får inte tolka singelmodellernas coef innan man har sett att de är sign. i multimodellen.
Nu när vi säkerställt att det är sign. på multinivå kan vi studera de univariata resultaten.
summary(fit1)## Response math :
##
## Call:
## lm(formula = math ~ dagis + specialdagis + tidextra + parents_educ +
## boy + kva1 + kva2 + kva3, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.6194 -3.9385 0.0385 3.5890 21.2494
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.30918 1.28959 18.850 < 2e-16 ***
## dagis 6.62620 1.04224 6.358 4.74e-10 ***
## specialdagis -0.42590 0.96831 -0.440 0.66025
## tidextra 1.74446 0.62757 2.780 0.00565 **
## parents_educ 0.10187 0.42559 0.239 0.81092
## boy -0.92928 0.56018 -1.659 0.09778 .
## kva1 5.55860 0.99201 5.603 3.53e-08 ***
## kva2 0.03919 0.92957 0.042 0.96639
## kva3 -0.97390 1.01739 -0.957 0.33892
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.044 on 484 degrees of freedom
## Multiple R-squared: 0.3278, Adjusted R-squared: 0.3167
## F-statistic: 29.5 on 8 and 484 DF, p-value: < 2.2e-16
##
##
## Response literacy :
##
## Call:
## lm(formula = literacy ~ dagis + specialdagis + tidextra + parents_educ +
## boy + kva1 + kva2 + kva3, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.5708 -3.9374 0.2228 4.0626 17.9856
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.6699 1.3185 11.884 < 2e-16 ***
## dagis 2.3705 1.0656 2.224 0.02658 *
## specialdagis -2.5767 0.9901 -2.603 0.00954 **
## tidextra 6.5216 0.6417 10.164 < 2e-16 ***
## parents_educ 2.1329 0.4351 4.901 1.3e-06 ***
## boy -1.4721 0.5728 -2.570 0.01046 *
## kva1 1.7917 1.0143 1.766 0.07795 .
## kva2 -0.7848 0.9504 -0.826 0.40934
## kva3 -0.6246 1.0402 -0.600 0.54846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.18 on 484 degrees of freedom
## Multiple R-squared: 0.3841, Adjusted R-squared: 0.3739
## F-statistic: 37.72 on 8 and 484 DF, p-value: < 2.2e-16
##
##
## Response social :
##
## Call:
## lm(formula = social ~ dagis + specialdagis + tidextra + parents_educ +
## boy + kva1 + kva2 + kva3, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.8041 -4.6777 0.7221 6.1959 18.5078
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.0798 1.7238 19.190 < 2e-16 ***
## dagis 8.7089 1.3932 6.251 8.95e-10 ***
## specialdagis -3.8211 1.2943 -2.952 0.00331 **
## tidextra 5.8011 0.8389 6.915 1.48e-11 ***
## parents_educ -0.3091 0.5689 -0.543 0.58713
## boy -1.8019 0.7488 -2.406 0.01649 *
## kva1 3.4355 1.3260 2.591 0.00986 **
## kva2 2.2143 1.2426 1.782 0.07536 .
## kva3 -0.9401 1.3599 -0.691 0.48974
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.079 on 484 degrees of freedom
## Multiple R-squared: 0.3054, Adjusted R-squared: 0.2939
## F-statistic: 26.6 on 8 and 484 DF, p-value: < 2.2e-16Ur summary(fit1) (längre ner i vår text) ser vi att R squared ligger på 0.3451, 0.3821, 0.3192, för math, literacy, social vilket är bra (högre än Ingers i föreläsningsexemplet).
Conclusions
Vi behöver ej titta på bakgrundsvariabler för att tolka. Vi kommer bara kommentera tolkningen på de signifikanta variablerna.
En sammanfattning av alla koef ges av
fit1##
## Call:
## lm(formula = cbind(math, literacy, social) ~ dagis + specialdagis +
## tidextra + parents_educ + boy + kva1 + kva2 + kva3, data = df)
##
## Coefficients:
## math literacy social
## (Intercept) 24.30918 15.66991 33.07975
## dagis 6.62620 2.37048 8.70888
## specialdagis -0.42590 -2.57666 -3.82115
## tidextra 1.74446 6.52164 5.80112
## parents_educ 0.10187 2.13287 -0.30911
## boy -0.92928 -1.47213 -1.80187
## kva1 5.55860 1.79167 3.43548
## kva2 0.03919 -0.78484 2.21431
## kva3 -0.97390 -0.62465 -0.94007De som spelar roll att tolka är de intressevariabler som är signifikanta. Vi plockar ut dessa ur tabellen ovan, och får följande tabell:
Coefficients:
math literacy social
(Intercept) 23.73265 16.07515 32.51889
dagis 7.40438 1.92834 9.49967
specialdagis EjSign. -2.41103 -3.19622
tidextra 1.52378 6.57021 5.55214Av den sistnämnda tabellen ser vi nu att vi kan göra tolkningar (vi studerade summary(fit1) för att se signifikans). Tolkningen är följande:
För literacy och social gäller att de som går i specialdagis, on average, har sämre testscore med -2.4 i literacy och -3.2 i social.
För alla test-scores gäller att
- De som går på dagis har du bättre test-score (hur mycket bättre ses av coef storlek: 7.4 för math 1.9 för literacy och 9.5 för social). Intressant är att minst ökning sker för literacy.
- De som spenderar över 15 timmar, i stället för 15 timmar, presterar on average bättre på provet (hur mycket bättre ses av coef storlek: 1.5 för math, 6.6 för literacy, 5.6 för social). Den koef som sticker ut är för math, där de som spenderar över 15 timmar bara har 1.5 bättre test-score on average, jämfört med koef för literacy och social där den ligger runt 6 dvs fyra gånger större.
Intercept tolkas sällan, det är den som tillåter oss att dra line of best fit, men här säger den oss något om vad test score är, on average, när samtliga variabler är satta till noll. Detta kan jämföras med histogrammen över våra Y, men är ej superviktig för tolkningen.
FRÅGA C
I val av metod ovan skrev vi följande:
Frågor att ställa oss innan vi bestämmer teknik (vi använder Ingers beslutsträd Enligt Lecture 1 slide 42):
- What type of relationships is being examined? Dependent or independent? Vi tycker dependent, för de tre Y vi har är sådant man typiskt sett vill se förklaras utav bakgrund. SVAR: DEPENDENT.
- Number of dependent. Vi har Y1 Y2 Y3 så SVAR: SEVERAL DEPENDENT IN SINGLE RELATIONSHIP.
- Measurement scale of the dependent variables. Våra Y går från 0 til 50. SVAR: NUMERICAL.
- Measurement scale of the explanatora variables. SVAR: Numerical -> cancorr och mutlivariate regression, SVAR: categorical -> manova. Vi har både numerical och categorical.
Då gruppstorlekarna är för små för att göra Manova (se kommentar tidigare också) samt att frågeställningen är “taking background variables into account” vilket regression hanterar mycket bra, så väljer vi att göra en multivariat multipel regression.
Om vi inte väljer denna metod, kan vi välja Canonical correlation eller Manova. Det vi är mest sugna på att använda är Manova, eftersom våra X är kategoriska - så vi skulle kunna analysera skillnader i medelvärden Y1 Y2 Y3 beroende på vilken grupp du är i (huruvida du är pojke/flicka, om du gick på dagis eller ej, samt om du gick på ett speciellt dagis elelr ej). Det är tre stycken ja/nej-frågor så vi skulle kunna skapa \(2^3=8\) grupper:
boy, dagis, specialdagis
100
110
101
111
000
010
001
011Eller så skulle vi kunna intressera oss för frågan “är pojkar duktigare än flickor?”, “spelar dagis roll eller ej för test-score??”, “inverkar specialdagis på test-score?” dvs testa en i taget.
Vi föredrar att använda regression eftersom vi då kan hålla bakgrundvariabler konstanta på ett intuitivt sätt (se ekvationen \(math+literacy+social = dagis+specialdagis+tidextra+parentseduc+boy+kva1+kva2+kva3\)) samt att vi kan tolka koef olika storlek på så sätt som vi gjorde under conclusion.
FRÅGA D
Frågan är 2) Is there any grouping of individuals that have certain combinations of variables in common, and if so: what is it that distinguishes these groups?
De klassificeringsmetoder vi har lärt oss i kursen är discriminant anlaysis och cluster analysis.
Vi väljer inte discriminant för att detta handlar om att classifiera in i grupper (t.ex. fisksort baserat på deras mått) och här har vi inte direkt någon klassificeringsvariabel som är intressant att titta på, givet frågan. Därmed är det cluster analysis som kan användas.
När används clustering? När vi vill gruppera objekt så att de inom gruppen är lika och de mellan grupper är olika. Exakt vad vi menar med lika beror på vilken typ av variabler vi är intresserade utav (för kontinuerliga har vi tre olika, för ordinal data ranking, för dummy proportion match). Ofta kan man se cluster med blotta ögat genom att göra en scatter.
Clustering passar vårt exempel eftersom frågan lyder “any grouping of individuals that have certain combinations of variables in common” dvs gruppera individer så att deras variabler liknar varandra.
Är det hierarkisk eller ej? Av formuleringen “any grouping” och “certain combinations” låter som att vi inte på förhand har en gissning om vilka som skulle höra ihop. Därmed är hierarkisk clustering det bästa valet.
Hur ska vi genomföra denna hierarkiska clustering? Hur vi ska göra den beror lite på vilken fråga vi vill svara på:
- Om vi vill se hur olika test scores hänger ihop (vi har ju sett att de har en korrelation och
pairs()ger intrycket av att ett samband finns) så kan vi gruppera ihop med avseende på Y. Givet att det finns en korerlation, och variablerna är kontinuerliga, kan vi använda distansmåtten euklides, mahalanobis, statistical dist. Då range av Y är ungefär samma och de är mätta på comparable scales så kan vi använda euklides avstånd (vilket är en fördel då standardisering kan påverka cluster lösningen). - Om vi vill se hur olika förklaringsvaribler (X) hänger ihop, kan vi klustra dem och läsa av dendogrammet.
- Om vi vill se hur olika intressevariabler hänger ihop kan vi göra som ovan. Då dessa är dummys är den lämpliga distance measure proportion matched.
- Om vi vill se hur olika bakgrundsvariabler hänger ihop skulle vi kanske vara intresserade av att clustra ihop olika parents education.
Dessvärre är inte clustering så bra på att fånga upp det som fråga (1) sökte, nämligen The relationship between knowledge variables and kindergarten experience, taking background variables into account. Om vi tror att Y beror av X så kommer clustering att udda resultat då den behandlar alla variabler som variabler - den gör ingen skillnad på X och Y. Därför har jag ovan exemplifierat på när man vill finna mönster mellan olika X (punkt 4, 3 och 2) eller mellan olika Y (punkt 1)
Detta är viktigt att ha i baktanke när vi tolkar variabeln \(tidextra\) som kommer att definieras snart - vi gör för att kunna göra en mer intuitiv tolkning en sak som kan lura den statistiskt okunnige då basline ej är den med flest observationer. Men då detta ej är grovt mycket fler, och läsaren till denna text är medveten om detta så är det inga större farhågor. Vi ville bara ge full disclusure så att tolkningen av coef i vår regression blir så som de ska vara.↩