<\body> <\input| >> x=rnorm(10000) <\input| >> hist(x,n=30);v() <\output> |ps>||||||> <\input| >> \; > How far are 5% of the points? <\input| >> x=sort(x) <\input| >> length(x)/100*5 <\output> [1] 500 <\input| >> x[500] <\output> [1] -1.627749 <\input| >> x[length(x)-500] <\output> [1] 1.629363 <\input| >> \; > We see that 5% of the points are 1.6 standard deviations below the mean, and 5% 1.6 above the mean. <\input| >> x[250] <\output> [1] -1.945368 <\input| >> x[length(x)-250] <\output> [1] 1.971970 <\input| >> \; > If we look at a distance of 1.96 std deviations away from the mean, or more, we'll see 5% of the points. There is a function that does this: <\input| >> qnorm(0.025) <\output> [1] -1.959964 <\input| >> qnorm(0.05) <\output> [1] -1.644854 <\input| >> qnorm(0.005) <\output> [1] -2.575829 <\input| >> \; > It is convenient to remember that if we are ~2 std. deviations away, we are at the 5% level, and at 3 std. deviations we are at less than 1% level. Even for this there is a function: <\input| >> pnorm(-3)*2 <\output> [1] 0.002699796 <\input| >> pnorm(-2)*2 <\output> [1] 0.04550026 <\input| >> pnorm(-1)*2 <\output> [1] 0.3173105 <\input| >> \; > Let us take 2 points that are taken from a normal, and look where they lie: <\input| >> x=rnorm(2000) <\input| >> y=rnorm(2000) <\input| >> plot(x,y,pch=".",xlim=c(-3,3),ylim=c(-3,3));v() <\output> |ps>||||||> <\input| >> grid(lwd=1,lty=1,col=2);v() <\output> |ps>||||||> <\input| >> \; > Let us count how many points are in each box <\input| >> x.cut=cut(x,-3:3) <\input| >> x.cut[1:20] <\output> \ [1] (0,1] \ \ (0,1] \ \ (1,2] \ \ (-2,-1] (0,1] \ \ (0,1] \ \ (1,2] \ \ (-1,0] \ (0,1] \ [10] (0,1] \ \ (-3,-2] (-1,0] \ (-2,-1] (-1,0] \ (-1,0] \ (0,1] \ \ (1,2] \ \ (-1,0]\ [19] (-1,0] \ (0,1] \ Levels: (-3,-2] (-2,-1] (-1,0] (0,1] (1,2] (2,3] <\input| >> \; > You see that each point was replaced by the interval that it is in. <\input| >> y.cut=cut(y,-3:3) <\input| >> tab.xy=table(x.cut,y.cut) <\input| >> tab.xy <\output> \ \ \ \ \ \ \ \ \ y.cut x.cut \ \ \ \ (-3,-2] (-2,-1] (-1,0] (0,1] (1,2] (2,3] \ \ (-3,-2] \ \ 1 \ \ \ \ \ \ 7 \ \ \ \ \ 15 \ \ \ \ 15 \ \ \ \ 9 \ \ \ \ 3 \ \ \ (-2,-1] \ \ 4 \ \ \ \ \ 31 \ \ \ \ \ 86 \ \ \ \ 85 \ \ \ 37 \ \ \ \ 8 \ \ \ (-1,0] \ \ 14 \ \ \ \ \ 82 \ \ \ \ 238 \ \ \ 242 \ \ \ 96 \ \ \ 13 \ \ \ (0,1] \ \ \ 16 \ \ \ \ \ 95 \ \ \ \ 226 \ \ \ 234 \ \ \ 83 \ \ \ 16 \ \ \ (1,2] \ \ \ \ 4 \ \ \ \ \ 34 \ \ \ \ \ 94 \ \ \ \ 99 \ \ \ 40 \ \ \ \ 7 \ \ \ (2,3] \ \ \ \ 3 \ \ \ \ \ \ 9 \ \ \ \ \ 13 \ \ \ \ 24 \ \ \ \ 4 \ \ \ \ 2 \ <\input| >> \; > Now we know how many are in each cell. Let us plot these: <\input| >> library(lattice) <\input| >> levelplot(tab.xy);v() <\output> |ps>||||||> <\input| >> \; > Now let us do the same, but with a finer grid, and more points: <\input| >> x=rnorm(1000000); y=rnorm(1000000) <\input| >> x.cut=cut(x, seq(-3,3, length=50) ); y.cut=cut(y, seq(-3,3,length=50) ) <\input| >> tab.xy=table(x.cut,y.cut) <\input| >> levelplot(tab.xy,colorkey=F,col.regions=rainbow(100));v() <\output> |ps>||||||> <\input| >> \; > Here we see why it makes sense to talk about sum of squares for normaly distributed data: All points for which +y=constant>, i.e. points on a circle, are equally likely. That is not true for other distributions. For example, the log-normal distribution: <\input| >> hist(exp(x),br=seq(0,1000,by=0.05),xlim=c(0,1.5));v() <\output> |ps>||||||> <\input| >> \; > Let us do a 2d plot in the same way here: <\input| >> x.cut=cut(exp(x), seq(0,1.5, length=50) ); y.cut=cut(exp(y), seq(0,1.5,length=50) ) <\input| >> tab.xy=table(x.cut,y.cut) <\input| >> levelplot(tab.xy,colorkey=F,col.regions=rainbow(100));v() <\output> |ps>||||||> <\input| >> \; > This is why it makes sense to try to minimize sum of squares for nomaly distributed data (with equal variances), but not always for non-normaly distributed data. <\input| >> a=read.table("data/regression.txt",head=T) <\input| >> plot( growth ~ tannin, data=a, pch=16);v() <\output> |ps>||||||> <\input| >> reg=lm( growth ~ tannin, data=a) <\input| >> summary(reg) <\output> \; Call: lm(formula = growth ~ tannin, data = a) \; Residuals: \ \ \ \ Min \ \ \ \ \ 1Q \ Median \ \ \ \ \ 3Q \ \ \ \ Max\ -2.4556 -0.8889 -0.2389 \ 0.9778 \ 2.8944\ \; Coefficients: \ \ \ \ \ \ \ \ \ \ \ \ Estimate Std. Error t value Pr(\\|t\|) \ \ \ (Intercept) \ 11.7556 \ \ \ \ 1.0408 \ 11.295 9.54e-06 *** tannin \ \ \ \ \ \ -1.2167 \ \ \ \ 0.2186 \ -5.565 0.000846 *** --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ \; Residual standard error: 1.693 on 7 degrees of freedom Multiple R-Squared: 0.8157, \ \ \ \ Adjusted R-squared: 0.7893\ F-statistic: 30.97 on 1 and 7 DF, \ p-value: 0.000846\ <\input| >> \; > Now let us look at aids to help us see how good the model is: <\input| >> two.by.two=t(matrix(1:4,2,2)) <\input| >> layout(two.by.two); par(cex=0.7) <\input| >> plot(reg);v(width=6,height=6) <\output> |ps>||||||> <\input| >> layout(1) <\input| >> \; > Let us look what these plots look like when a linear fit is not a good approximation: <\input| >> x=1:50 <\input| >> y=x^2+rnorm(50,sd=100) <\input| >> plot( y ~ x );v() <\output> |ps>||||||> <\input| >> reg=lm( y ~ x ) <\input| >> abline(reg,col=2); v() <\output> |ps>||||||> <\input| >> summary(reg) <\output> \; Call: lm(formula = y ~ x) \; Residuals: \ \ \ \ Min \ \ \ \ \ 1Q \ Median \ \ \ \ \ 3Q \ \ \ \ Max\ -335.50 -191.21 \ -49.87 \ 137.18 \ 569.01\ \; Coefficients: \ \ \ \ \ \ \ \ \ \ \ \ Estimate Std. Error t value Pr(\\|t\|) \ \ \ (Intercept) -484.816 \ \ \ \ 69.086 \ -7.018 6.91e-09 *** x \ \ \ \ \ \ \ \ \ \ \ \ 51.827 \ \ \ \ \ 2.358 \ 21.981 \ \ 2e-16 *** --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ \; Residual standard error: 240.6 on 48 degrees of freedom Multiple R-Squared: 0.9096, \ \ \ \ Adjusted R-squared: 0.9077\ F-statistic: 483.1 on 1 and 48 DF, \ p-value: \ 2.2e-16\ <\input| >> layout(two.by.two) <\input| >> plot(reg); v(width=6,height=6) <\output> |ps>||||||> <\input| >> layout(1) <\input| >> \; > Another example: <\input| >> x=1:50 <\input| >> y=x+rnorm(50,sd=1:50/3) <\input| >> plot(y~x);v() <\output> |ps>||||||> <\input| >> reg=lm(y~x) <\input| >> abline(reg,col=2);v() <\output> |ps>||||||> <\input| >> layout(two.by.two); plot(reg); v( width=6,height=6); layout(1) <\output> |ps>||||||> <\input| >> \; > Let us do the following: <\input| >> x=1:2 <\input| >> y=rnorm(2) <\input| >> plot(x,y);v() <\output> |ps>||||||> <\input| >> reg=lm(y~x) <\input| >> abline(reg,col=2);v() <\output> |ps>||||||> <\input| >> \; > Oh, nice! Perfect fit! Compare this with the fit we would get if we had a slope of 0: <\input| >> reg0=lm(y~1) <\input| >> reg0 <\output> \; Call: lm(formula = y ~ 1) \; Coefficients: (Intercept) \ \ \ \ \ 0.03782 \ <\input| >> abline(0.03782,0,col=3);v() <\output> |ps>||||||> <\input| >> \; > Isn't the red line a much better fit than the green one? No, it isn't! We can fit a line through any two points! So, if we have two measurements, it does not make sense to fit more than one parameter. <\input| >> anova(reg,reg0) <\output> Analysis of Variance Table \; Model 1: y ~ x Model 2: y ~ 1 \ \ Res.Df \ \ \ \ RSS Df Sum of Sq F Pr(\F) 1 \ \ \ \ \ 0 \ 0.0000 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ 1 \ 1.6424 -1 \ \ -1.6424 \ \ \ \ \ \ \ \ <\input| >> \; > These are the degrees of freedom of the model. Every independent measurement adds a degree of freedom, and every (independent) fitted parameter takes one degree away. Let us do this again <\input| >> x=1:3 <\input| >> y=rnorm(3) <\input| >> plot(x,y);v() <\output> |ps>||||||> <\input| >> reg=lm(y~x);abline(reg,col=2);v() <\output> |ps>||||||> <\input| >> reg0=lm(y~1);abline(reg0$coef[1],0,col=3);v() <\output> |ps>||||||> <\input| >> \; > So, again - the red line looks like a much better fit than the green one - but is it? Was it worth it to add another coeficient? Is the added fit significant? <\input| >> anova(reg,reg0) <\output> Analysis of Variance Table \; Model 1: y ~ x Model 2: y ~ 1 \ \ Res.Df \ \ \ \ \ RSS Df Sum of Sq \ \ \ \ \ F Pr(\F) 1 \ \ \ \ \ 1 \ 0.56735 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ 2 \ 1.55077 -1 \ -0.98342 1.7333 0.4135 <\input| >> \; > Analysis of variance tells us here than adding another coefficient did not make the model better, because the p-value is pretty big. Let us look at an example (from the book...) <\input| >> a=read.table("data/factorial.txt",head=T) <\input| >> a <\output> \ \ \ growth diet \ coat 1 \ \ \ \ 6.6 \ \ \ A light 2 \ \ \ \ 7.2 \ \ \ A light 3 \ \ \ \ 6.9 \ \ \ B light 4 \ \ \ \ 8.3 \ \ \ B light 5 \ \ \ \ 7.9 \ \ \ C light 6 \ \ \ \ 9.2 \ \ \ C light 7 \ \ \ \ 8.3 \ \ \ A \ dark 8 \ \ \ \ 8.7 \ \ \ A \ dark 9 \ \ \ \ 8.1 \ \ \ B \ dark 10 \ \ \ 8.5 \ \ \ B \ dark 11 \ \ \ 9.1 \ \ \ C \ dark 12 \ \ \ 9.0 \ \ \ C \ dark <\input| >> layout(matrix(1:2,1,2));plot( growth ~ diet, data=a); <\input| >> plot( growth ~ coat, data=a);v();layout(1) <\output> |ps>||||||> > Now we can try to fit linear models: <\input| >> reg1=lm(growth~coat,data=a) <\input| >> reg2=lm(growth~diet+coat,data=a) <\input| >> reg3=lm(growth~coat:diet+coat+diet,data=a) <\input| >> anova(reg1,reg2,reg3) <\output> Analysis of Variance Table \; Model 1: growth ~ coat Model 2: growth ~ diet + coat Model 3: growth ~ coat:diet + coat + diet \ \ Res.Df \ \ \ RSS Df Sum of Sq \ \ \ \ \ F \ Pr(\F) \ 1 \ \ \ \ 10 5.5167 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ 8 2.8567 \ 2 \ \ \ 2.6600 3.6774 0.09069 . 3 \ \ \ \ \ 6 2.1700 \ 2 \ \ \ 0.6867 0.9493 0.43833 \ --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > It is easier to call the function aov, which will call lm along the formula: <\input| >> model1=aov( growth ~ coat+diet+diet:coat, data=a) > diet*coat is a shorthand for coat+diet+diet:coat \; <\input| >> model1=aov(growth~coat*diet,data=a) <\input| >> summary(model1) <\output> \ \ \ \ \ \ \ \ \ \ \ \ Df \ Sum Sq Mean Sq F value \ Pr(\F) \ coat \ \ \ \ \ \ \ \ 1 2.61333 2.61333 \ 7.2258 0.03614 * diet \ \ \ \ \ \ \ \ 2 2.66000 1.33000 \ 3.6774 0.09069 . coat:diet \ \ \ 2 0.68667 0.34333 \ 0.9493 0.43833 \ Residuals \ \ \ 6 2.17000 0.36167 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> layout(two.by.two);par(cex=0.8);plot(model1) <\input| >> v(width=6,height=6);layout(1) <\output> |ps>||||||> <\input| >> model2=update(model1,~.-diet:coat) <\input| >> anova(model1,model2) <\output> Analysis of Variance Table \; Model 1: growth ~ coat + diet + diet:coat Model 2: growth ~ coat + diet \ \ Res.Df \ \ \ \ \ RSS Df Sum of Sq \ \ \ \ \ F Pr(\F) 1 \ \ \ \ \ 6 \ 2.17000 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ 8 \ 2.85667 -2 \ -0.68667 0.9493 0.4383 <\input| >> \; > So, it says that model1 is not significantly better than model2. <\input| >> summary(model2) <\output> \ \ \ \ \ \ \ \ \ \ \ \ Df \ Sum Sq Mean Sq F value \ Pr(\F) \ coat \ \ \ \ \ \ \ \ 1 2.61333 2.61333 \ 7.3186 0.02685 * diet \ \ \ \ \ \ \ \ 2 2.66000 1.33000 \ 3.7246 0.07190 . Residuals \ \ \ 8 2.85667 0.35708 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > Let us remove diet: <\input| >> model3=update(model2,~.-diet) <\input| >> anova(model2,model3) <\output> Analysis of Variance Table \; Model 1: growth ~ coat + diet Model 2: growth ~ coat \ \ Res.Df \ \ \ \ RSS Df Sum of Sq \ \ \ \ \ F Pr(\F) \ 1 \ \ \ \ \ 8 \ 2.8567 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ 10 \ 5.5167 -2 \ \ -2.6600 3.7246 0.0719 . --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > So, model2 is not significantly better than model3. <\input| >> summary(model3) <\output> \ \ \ \ \ \ \ \ \ \ \ \ Df Sum Sq Mean Sq F value \ Pr(\F) \ coat \ \ \ \ \ \ \ \ 1 2.6133 \ 2.6133 \ 4.7372 0.05457 . Residuals \ \ 10 5.5167 \ 0.5517 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > Now coat is not significant either! <\input| >> layout(matrix(1:2,1,2));plot( growth ~ diet, data=a); <\input| >> plot( growth ~ coat, data=a);v();layout(1) <\output> |ps>||||||> <\input| >> \; > We see that diet C seems very different from diet A and B, but that they are not very different from each other. Let us just make a model that depends on whether the diet is C or not. <\input| >> dietC=as.factor(a$diet=="C") <\input| >> model4=update(model3,~.+dietC) <\input| >> anova(model3,model4) <\output> Analysis of Variance Table \; Model 1: growth ~ coat Model 2: growth ~ coat + dietC \ \ Res.Df \ \ \ RSS Df Sum of Sq \ \ \ \ \ F \ Pr(\F) \ 1 \ \ \ \ 10 5.5167 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ 9 2.9817 \ 1 \ \ \ 2.5350 7.6518 0.02189 * --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > So, model4 is significantly better than model3. Let us try adding interaction: <\input| >> model5=update(model4,~.+dietC:coat) <\input| >> anova(model4,model5) <\output> Analysis of Variance Table \; Model 1: growth ~ coat + dietC Model 2: growth ~ coat + dietC + coat:dietC \ \ Res.Df \ \ \ \ RSS Df Sum of Sq \ \ \ \ \ F Pr(\F) 1 \ \ \ \ \ 9 2.98167 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ 8 2.70000 \ 1 \ \ 0.28167 0.8346 0.3877 <\input| >> \; > This addition is not worthwhile. \; <\initial> <\collection> <\references> <\collection> |?>> |?>> |?>> |?>> |?>> <\auxiliary> <\collection> <\associate|toc> Lecture 6 - Linear regression and analysis of variance |Some more about the normal distribution> |2 dimensional distributions>