Aanwijzingen bij het maken van de sommen in R uit het boek van Rice ------------------------------------------- Algemeen ------------------------------------------- Om een bepaalde functie op de rijen of kolommen van een matrix toe te passen kun je de functie apply gebruiken: > apply(x,2,mean) en > apply(x,1,mean) geven respectievelijk de kolomsommen en de rijsommen van matrix x. Een selectie van bepaalde kolommen of rijen van een matrix x kan als volgt worden gemaakt: > x[2:5,] geeft rij 2 t/m 5 van matrix x. > x[,c(1,4,5,6)] geeft kolom 1,4,5 en 6 van matrix x. ------------------------------------------- Aanwijzingen bij Hoofdstuk 14 ------------------------------------------- Enkelvoudige lineaire regressie "y_i = beta_0 + beta_1 * x_i + e_i" kun je uitrekenen met de functie > z=lsfit(x,y) waar x en y vectoren zijn. Om zonder beta_0 alleen "y_i = beta_1 * x_i + e_i" te fitten gebruik je > z=lsfit(x,y,intercept=F) In z staat allerlei informatie, bijvoorbeeld > z$coef geeft de geschatte parameters. > z$res geeft de geschatte residuen. In ls.diag(z) vind je allerlei informatie omtrent de schatting, bijvoorbeeld > ls.diag(z)$std.dev geeft schatting van sigma > ls.diag(z)$std.err geeft de geschatte standaardafwijkingen van de parameter- schattingen. Meervoudige lineaire regressie kan met dezelfde functie lsfit, maar nu moet als tweede argument de design matrix worden toegevoegd. Voor meervoudige lineaire regressie "y_i = beta_0 + beta_1 * x_[i,1] + ...+ beta_[p-1] * x_[i,p-1] + e_i" moet je eerst de design matrix aanmaken. Als je de kolommen van de design matrix hebt opgeslagen in x1,x2,...,xp dan maak je de matrix met > X=cbind(x1,x2,...,xp) of met > X=matrix(c(x1,x2,...,xp),ncol=p,nrow=n) Vervolgens kun je weer lsfit gebruiken: > z=lsfit(X,y) waar nu X de design matrix is en y de waarnemingen van de afhankelijke variabele. z$coef, z$res, ls.diag(z)$std.dev, ls.diag(z)$std.err werken analoog als bij enkelvoudige lineaire regressie. Om scatterplots te maken van de waarnemingen van alle variabelen (onafhankelijke en afhankelijke beide) kun je de functie pairs(matrix) gebruiken. Wil je slechts een selectie van de scatterplots dan kun je bijvoorbeeld de tweede tegen de vierde kolom van X plotten met > plot(X[,2],X[,4]). ------------------------------------------- Aanwijzingen bij Hoofdstuk 12 ------------------------------------------- Variantie analyse werkt in R met de functie aov(). Om deze functie aan te roepen heb je een "data.frame" nodig. Om je data.frame te maken zet je alle metingen (dus de y-waarden van alle groepen van alle factoren) eerst in een vector (totaledata) > totaledata=c(.......) Vervolgens maak je per factor een vector "factor" aan die aan alle elementen uit de vector totaledata de naam meegeeft van de groep van die factor. De volgorde van de data in "totaledata" bepaalt de volgorde van de elementen in de "factor" vectoren. Bijvoorbeeld als je een medicijnfactor hebt, waarbij je vier soorten medicijnen hebt, A,B,C en D, dan kan je factor "medicijn" er als volgt uitzien: > medicijn=factor(c(rep(1,10),rep(2,10),rep(3,10),rep(4,10)),labels=c("A","B","C","D")) als je per groep 10 patienten hebt gemeten. Let op dat je de juiste volgorde neemt, het kan ook zijn dat je > medicijn=factor(rep(1:4,10),labels=c("A","B","C","D")) nodig hebt, dat hangt van de indeling in rijen en kolommen af. Dit kun je controleren door medicijn en totaledata beide op te vragen en te kijken of het matcht. Je kunt meerdere factoren aanmaken op deze manier. Vervolgens maak je het data.frame met (voorbeeld met 2 factoren): > mijn.data = data.frame(totaledata,factor1,factor2) De variantie-analyse gaat dan als volgt voor 1-factor model: > analyse=aov(totaledata~factor1,mijn.data) en voor 2-factor model > analyse=aov(totaledata~factor1+factor2,mijn.data) voor 2-factor model met interacties > analyse=aov(totaledata~factor1+factor2+factor1:factor2,mijn.data) of kort: > analyse=aov(totaledata~factor1*factor2,mijn.data) voor 2-factor model met alleen factor 1 en interacties > analyse=aov(totaledata~factor1+factor1:factor2,mijn.data) enzovoort. Met > summary(analyse) krijg je een tabel als in Example A op p.483 in Rice Met > model.tables(analyse) krijg je de geschatte alfa_i (die sommeren tot 0) en met > model.tables(analyse,type="means") krijg je de geschatte (mu+alfa_i) voor alle groepen (die sommeren tot mu) en ook de geschatte waarde voor mu in "grand mean". In de methode van Tukey heb je kwantielen van de studentized range distribtution nodig. Die vind je in R met > qtukey(alfa,df1,df2) De Kruskal-Wallis toets werkt alleen voor het 1-factor model. > kruskal.test(totaledata,factor1) geeft in de output de waarde van de chikwadraat-statistiek en de overschrijdingskans (p-value). ------------------------------------------- Aanwijzingen bij Hoofdstuk 13 ------------------------------------------- Om Fisher's exacte toets uit te voeren, gebruik je fisher.test() in R. Om de chikwadraat toets uit te voeren, gebruik je chisq.test() in R. De functie chisq.test gebruikt standaard de Yates correctie voor 2x2 tabellen. Dit moet je uitzetten; je kunt dit doen door het argument correct=F mee te geven. Controleer in geval van een experiment met weinig waarnemingen of je voldoet aan de conditie van de toepasbaarheid van de chikwadraattoets: onder H0 moet voor elke cel de verwachting van het aantal waarnemingen in de cel minstens 5 bedragen.