HCO16, dealing with heterogeneity and error bars
General linear model, tot nu toe hebben we steeds general linear models gebruikt. Deze maken
gebruik van ordinary least squares (OLS) methodes. Voorbeelden zijn regressie, ANOVA en ANCOVA.
Allemaal kunnen ze uitgevoerd worden middels de lm() functie. Er is echter een probleem, want deze
methodes zijn beperkt tot:
- Data met gelijke variantie in de residuen, als dat niet geval is en je variantie dus heterogeen
is, kan je ervoor kiezen om gls (generalized least squares) te gebruiken.
- Data met normaal verdeelde residuen, als dat niet het geval is, kan je een generalized linear
model (glm) gebruiken.
Data, we zullen gls gaan behandelen aan de hand van
een voorbeeld waarvoor data van BioClive gebruikt
wordt. Er zijn meerdere plotjes die allemaal in 4
subplots onderverdeeld kunnen worden, zoals rechts
te zien is. In totaal zijn er 12 monoculturen, 12
plotten met 4 plantensoorten, 12 plotten met 8
plantensoorten, 8 plotten met 12 plantensoorten. In totaal heb je dus (3*12+8=)44 plotjes. Verder
zijn er 2 blocks dus daarmee kom je uit op 88 plotten in totaal en dus (88*4=)352 subplotjes. De
biodiversiteit van deze plotjes wordt gemanipuleerd door
ongewilde planten weg te halen, aangezien in dit
onderzoek gekeken wordt wat er met een ecosysteem
gebeurt als de biodiversiteit afneemt. Een paar vragen die
met dit experiment beantwoord kunnen worden, zijn:
- Biedt hogere plantdiversiteit een hogere weerstand
tegen onkruid invasie? Dit is te meten door te
kijken hoeveel onkruid per diversiteitslevels weg is
gehaald. De niche hypothese stelt dat meer soorten
samen meer niches invullen en invasie dan dus
minder makkelijk op kan treden. Je ziet rechts in de
bovenste grafiek dat dit inderdaad lijkt te kloppen.
- Is deze weerstand tegen onkruid invasie afhankelijk
van biomassa? Hoe hoger de biomassa is, hoe meer
fysieke ruimte al ingenomen is en hoe meer
lichtcompetitie er op zal treden met het onkruid.
- Is de relatie tussen gezaaide plant biomassa en
onkruid invasie afhankelijk van de plantdiversiteit?
Structuur data, rechts zie in een tabel een deel van de data
aangegeven en we hebben dus de volgende variabelen:
- Weed_weight, is hier de continue respons variabele (y)
- AGbiomass_weight, is een continue explanatory variabele (x)
- Div, is een categoriale explanatory variabele, factor met 4 levels (X)
Je hebt hier dus data om een ANCOVA mee uit te voeren. Voor dit voorbeeld
zijn we alleen geïnteresseerd in de main effects, omdat dat het iets
makkelijker maakt en ons model ziet er dus als
volgt uit: lm(y~x + X). Rechts zie je een
scatter plot die als het ware de ANCOVA
weergeeft: voor elke groep een regressie. Er is
echter een probleem, want er worden
negatieve waardes voorspelt. Gelukkig is dit
gemakkelijk op te lossen door de log te
nemen.
, lm(), als je een ANCOVA uit wil voeren zal je waarschijnlijk
eerst lm() uitproberen. Rechts is te zien dat de
diagnostische plot er voor de BioClive niet heel goed uitziet.
De relatie is vrij lineair en de data is op zich wel normaal
verdeeld, maar er is overduidelijk sprake van
heteroscedasticity. Je kan dan nog best een summary en
anova tabel maken, maar deze output kan je niet
vertrouwen, aangezien niet aan de aannames voldaan
wordt.
Heteroscedasticity, je kan ervoor kiezen om allerlei
transformaties uit te voeren op de data, zodat wel aan de
aannames van lm() voldaan wordt, maar een andere optie
is gls.
gls(), je ziet links dat de gls() en lm()
plot exact hetzelfde zijn op het feit
na dat de gls() diagnostische plot
geen rode lijn bevat. Als je de anova
en summary tabel uitprint, zal je
zien dat deze ook vrijwel hetzelfde
zijn. De opmaak is misschien iets
anders, maar de berekende
waardes zijn hetzelfde. Gls() doet essentieel dus hetzelfde als lm(): ze fitten een lineair model met
een gaussian (normaal) verdeling. We weten echter dat we deze output niet kunnen vertrouwen,
want we hebben met heteroscedasticity te maken.
Homoscedasticity, het probleem is dus dat de variantie van de residuen toeneemt met een toename
in biomassa. Daarnaast verschilt de variantie van de residuen ook nog per level van diversiteit. lm
gaat echter uit van homoscedasticity. Een oplossing om dan toch met de BioClive data te kunnen
werken, is een model te gebruiken dat toestaat dat de variantie verandert als functie van de
explanatory variabelen. Een voorbeeld van zo’n model is gls.
Gls (generalized least squares), werkt eigenlijk volgens dezelfde methode als OLS, maar is minder
streng voor de aanname dat de variantie homogeen verdeeld moet zijn. Gls staat namelijk
heteroscedasticity toe en/of een correlatie van de variantie met een explanatory variabele. Dit doet
gls door het verschil in variantie te meten i.p.v. aan te nemen dat de variantie hetzelfde is. Gls
modelleert de verandering in variantie als een functie van een explanatory variabele. Dit doe je door
het ‘weights’ argument toe te voegen aan de gls() functie.
varIdent, als je weet dat je met ongelijke
variantie te maken hebt voor
verschillende levels van je categoriale
explanatory variabele, maakt varIdent
het mogelijk dat elk level van diversiteit
zijn eigen variantie heeft. Hiermee geef
je aan dat de variantie voor elk level
geschat wordt i.p.v. aangenomen wordt.
Rechts zie je hoe je varIdent toepast in ons voorbeeld waarbij diversiteit de categoriale explanatory
variabele is:
weights = varIdent(form = ~1|div)
Volgens meneer Hautier ziet de variantie er nu al beter uit, maar ik zie eerlijk gezegd geen verschil.
General linear model, tot nu toe hebben we steeds general linear models gebruikt. Deze maken
gebruik van ordinary least squares (OLS) methodes. Voorbeelden zijn regressie, ANOVA en ANCOVA.
Allemaal kunnen ze uitgevoerd worden middels de lm() functie. Er is echter een probleem, want deze
methodes zijn beperkt tot:
- Data met gelijke variantie in de residuen, als dat niet geval is en je variantie dus heterogeen
is, kan je ervoor kiezen om gls (generalized least squares) te gebruiken.
- Data met normaal verdeelde residuen, als dat niet het geval is, kan je een generalized linear
model (glm) gebruiken.
Data, we zullen gls gaan behandelen aan de hand van
een voorbeeld waarvoor data van BioClive gebruikt
wordt. Er zijn meerdere plotjes die allemaal in 4
subplots onderverdeeld kunnen worden, zoals rechts
te zien is. In totaal zijn er 12 monoculturen, 12
plotten met 4 plantensoorten, 12 plotten met 8
plantensoorten, 8 plotten met 12 plantensoorten. In totaal heb je dus (3*12+8=)44 plotjes. Verder
zijn er 2 blocks dus daarmee kom je uit op 88 plotten in totaal en dus (88*4=)352 subplotjes. De
biodiversiteit van deze plotjes wordt gemanipuleerd door
ongewilde planten weg te halen, aangezien in dit
onderzoek gekeken wordt wat er met een ecosysteem
gebeurt als de biodiversiteit afneemt. Een paar vragen die
met dit experiment beantwoord kunnen worden, zijn:
- Biedt hogere plantdiversiteit een hogere weerstand
tegen onkruid invasie? Dit is te meten door te
kijken hoeveel onkruid per diversiteitslevels weg is
gehaald. De niche hypothese stelt dat meer soorten
samen meer niches invullen en invasie dan dus
minder makkelijk op kan treden. Je ziet rechts in de
bovenste grafiek dat dit inderdaad lijkt te kloppen.
- Is deze weerstand tegen onkruid invasie afhankelijk
van biomassa? Hoe hoger de biomassa is, hoe meer
fysieke ruimte al ingenomen is en hoe meer
lichtcompetitie er op zal treden met het onkruid.
- Is de relatie tussen gezaaide plant biomassa en
onkruid invasie afhankelijk van de plantdiversiteit?
Structuur data, rechts zie in een tabel een deel van de data
aangegeven en we hebben dus de volgende variabelen:
- Weed_weight, is hier de continue respons variabele (y)
- AGbiomass_weight, is een continue explanatory variabele (x)
- Div, is een categoriale explanatory variabele, factor met 4 levels (X)
Je hebt hier dus data om een ANCOVA mee uit te voeren. Voor dit voorbeeld
zijn we alleen geïnteresseerd in de main effects, omdat dat het iets
makkelijker maakt en ons model ziet er dus als
volgt uit: lm(y~x + X). Rechts zie je een
scatter plot die als het ware de ANCOVA
weergeeft: voor elke groep een regressie. Er is
echter een probleem, want er worden
negatieve waardes voorspelt. Gelukkig is dit
gemakkelijk op te lossen door de log te
nemen.
, lm(), als je een ANCOVA uit wil voeren zal je waarschijnlijk
eerst lm() uitproberen. Rechts is te zien dat de
diagnostische plot er voor de BioClive niet heel goed uitziet.
De relatie is vrij lineair en de data is op zich wel normaal
verdeeld, maar er is overduidelijk sprake van
heteroscedasticity. Je kan dan nog best een summary en
anova tabel maken, maar deze output kan je niet
vertrouwen, aangezien niet aan de aannames voldaan
wordt.
Heteroscedasticity, je kan ervoor kiezen om allerlei
transformaties uit te voeren op de data, zodat wel aan de
aannames van lm() voldaan wordt, maar een andere optie
is gls.
gls(), je ziet links dat de gls() en lm()
plot exact hetzelfde zijn op het feit
na dat de gls() diagnostische plot
geen rode lijn bevat. Als je de anova
en summary tabel uitprint, zal je
zien dat deze ook vrijwel hetzelfde
zijn. De opmaak is misschien iets
anders, maar de berekende
waardes zijn hetzelfde. Gls() doet essentieel dus hetzelfde als lm(): ze fitten een lineair model met
een gaussian (normaal) verdeling. We weten echter dat we deze output niet kunnen vertrouwen,
want we hebben met heteroscedasticity te maken.
Homoscedasticity, het probleem is dus dat de variantie van de residuen toeneemt met een toename
in biomassa. Daarnaast verschilt de variantie van de residuen ook nog per level van diversiteit. lm
gaat echter uit van homoscedasticity. Een oplossing om dan toch met de BioClive data te kunnen
werken, is een model te gebruiken dat toestaat dat de variantie verandert als functie van de
explanatory variabelen. Een voorbeeld van zo’n model is gls.
Gls (generalized least squares), werkt eigenlijk volgens dezelfde methode als OLS, maar is minder
streng voor de aanname dat de variantie homogeen verdeeld moet zijn. Gls staat namelijk
heteroscedasticity toe en/of een correlatie van de variantie met een explanatory variabele. Dit doet
gls door het verschil in variantie te meten i.p.v. aan te nemen dat de variantie hetzelfde is. Gls
modelleert de verandering in variantie als een functie van een explanatory variabele. Dit doe je door
het ‘weights’ argument toe te voegen aan de gls() functie.
varIdent, als je weet dat je met ongelijke
variantie te maken hebt voor
verschillende levels van je categoriale
explanatory variabele, maakt varIdent
het mogelijk dat elk level van diversiteit
zijn eigen variantie heeft. Hiermee geef
je aan dat de variantie voor elk level
geschat wordt i.p.v. aangenomen wordt.
Rechts zie je hoe je varIdent toepast in ons voorbeeld waarbij diversiteit de categoriale explanatory
variabele is:
weights = varIdent(form = ~1|div)
Volgens meneer Hautier ziet de variantie er nu al beter uit, maar ik zie eerlijk gezegd geen verschil.
