<\body> \ 1a. Use the data in the file factorial.txt. Compare the mean of growth when the coat is light to the mean when the coat is dark. Is the variance in the two cases significantly different? Are the distributions normal? <\input|r] > \; <\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| >> i=a$coat=="dark" <\input| >> var.test(a$growth[i],a$growth[!i]) <\output> \; \ \ \ \ \ \ \ \ F test to compare two variances \; data: \ a$growth[i] and a$growth[!i]\ F = 0.1618, num df = 5, denom df = 5, p-value = 0.06724 alternative hypothesis: true ratio of variances is not equal to 1\ 95 percent confidence interval: \ 0.02264239 1.15636435\ sample estimates: ratio of variances\ \ \ \ \ \ \ \ \ \ 0.1618112\ <\input| >> t.test(a$growth[i],a$growth[!i],var.eq=T) <\output> \; \ \ \ \ \ \ \ \ Two Sample t-test \; data: \ a$growth[i] and a$growth[!i]\ t = 2.1765, df = 10, p-value = 0.05457 alternative hypothesis: true difference in means is not equal to 0\ 95 percent confidence interval: \ -0.02214314 \ 1.88880980\ sample estimates: mean of x mean of y\ \ 8.616667 \ 7.683333\ <\input| >> b=a$growth[i] <\input| >> b=b-mean(b) <\input| >> b=b/sd(b) <\input| >> ks.test(b,pnorm) <\output> \; \ \ \ \ \ \ \ \ One-sample Kolmogorov-Smirnov test \; data: \ b\ D = 0.1693, p-value = 0.9954 alternative hypothesis: two.sided\ <\input| >> b=a$growth[!i];b=b-mean(b);b=b/sd(b) <\input| >> ks.test(b,pnorm) <\output> \; \ \ \ \ \ \ \ \ One-sample Kolmogorov-Smirnov test \; data: \ b\ D = 0.19, p-value = 0.9819 alternative hypothesis: two.sided\ <\input| >> \; > 1b. Do a linear model of growth, that takes into account coat color. Is it significantly better than not taking coat color into account? Compare to the p-value you got in 1a. <\input| >> reg=lm(growth~coat,data=a) <\input| >> reg <\output> \; Call: lm(formula = growth ~ coat, data = a) \; Coefficients: (Intercept) \ \ \ coatlight \ \ \ \ \ \ 8.6167 \ \ \ \ \ -0.9333 \ <\input| >> anova(reg) <\output> Analysis of Variance Table \; Response: growth \ \ \ \ \ \ \ \ \ \ 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| >> reg1=lm(growth~1,data=a) <\input| >> anova(reg,reg1) <\output> Analysis of Variance Table \; Model 1: growth ~ coat Model 2: growth ~ 1 \ \ Res.Df \ \ \ \ RSS Df Sum of Sq \ \ \ \ \ F \ Pr(\F) \ 1 \ \ \ \ 10 \ 5.5167 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ 11 \ 8.1300 -1 \ \ -2.6133 4.7372 0.05457 . --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > Using coat is borderline insignificant at a 5% level, and we get exactly the same p-value as above. \; 2a. Draw 20 random points from a normal, and another 2 from a normal with mean 0.5. Does a t-test see them as significantly different? Does a Wilcoxon test see them as significantly different? <\input| >> x=rnorm(20);y=rnorm(20,mean=0.5) <\input| >> t.test(x,y,var.eq=T) <\output> \; \ \ \ \ \ \ \ \ Two Sample t-test \; data: \ x and y\ t = -1.1355, df = 38, p-value = 0.2633 alternative hypothesis: true difference in means is not equal to 0\ 95 percent confidence interval: \ -1.0247315 \ 0.2882662\ sample estimates: mean of x mean of y\ 0.1317953 0.5000279\ <\input| >> wilcox.test(x,y) <\output> \; \ \ \ \ \ \ \ \ Wilcoxon rank sum test \; data: \ x and y\ W = 153, p-value = 0.2110 alternative hypothesis: true mu is not equal to 0\ <\input| >> \; > 2b. Do the same as you did in 2a 1000 times with the t-test, and 1000 times with the Wilcoxon test. Are the success rates of the two tests significantly different from each other? try to find at what distance (i.e. mean of the second 20 samples) they become different. <\input| >> compare=function(mean){\ \ x=matrix(rnorm(1000*20),1000,20) \ y=matrix(rnorm(1000*20,mean=mean),1000,20) tp=apply(cbind(x,y),1,function(xy) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t.test(xy[1:20],xy[21:40],var.eq=T)$p.value) \ wp=apply(cbind(x,y),1,function(xy) wilcox.test(xy[1:20],xy[21:40])$p.value) c(sum(tp\0.05),sum(wp\0.05)) } <\input| >> \; <\input| >> prop.test(compare(0.1),n=c(1000,1000)) <\output> \; \ \ \ \ \ \ \ \ 2-sample test for equality of proportions without continuity \ \ \ \ \ \ \ \ correction \; data: \ compare(0.1) out of c(1000, 1000)\ X-squared = 0, df = 1, p-value = 1 alternative hypothesis: two.sided\ 95 percent confidence interval: \ -0.02191499 \ 0.02191499\ sample estimates: prop 1 prop 2\ \ 0.067 \ 0.067\ <\input| >> x=sapply(seq(0.01,3,length=20),compare) <\input| >> x <\output> \ \ \ \ \ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [1,] \ \ 50 \ \ 89 \ 177 \ 296 \ 495 \ 688 \ 853 \ 918 \ 975 \ \ 995 \ \ 997 \ 1000 \ 1000 \ 1000 [2,] \ \ 52 \ \ 87 \ 163 \ 284 \ 472 \ 661 \ 825 \ 912 \ 964 \ \ 989 \ \ 997 \ 1000 \ 1000 \ 1000 \ \ \ \ \ [,15] [,16] [,17] [,18] [,19] [,20] [1,] \ 1000 \ 1000 \ 1000 \ 1000 \ 1000 \ 1000 [2,] \ 1000 \ 1000 \ 1000 \ 1000 \ 1000 \ 1000 <\input| >> plot(x[1,]);points(x[2,],col=2);v() <\output> |ps>||||||> <\input| >> \; > 3. Read the file strange.txt. Plot z vs. x. Do a linear regression of z vs. x. Does z increase or dercrease with x? Now do an analysis of variance of linear models for z. In the best model, does z increase or decrease with x? <\input| >> a=read.table("strange.txt",head=T,sep=",") <\input| >> reg1=lm(z~x,data=a) <\input| >> plot(z~x,data=a) <\input| >> abline(reg1);v() <\output> |ps>||||||> <\input| >> \; > There seems to be a positive slope <\input| >> reg1 <\output> \; Call: lm(formula = z ~ x, data = a) \; Coefficients: (Intercept) \ \ \ \ \ \ \ \ \ \ \ x \ \ \ \ \ \ 0.8757 \ \ \ \ \ \ 1.0115 \ <\input| >> reg2=aov(z~x*y,data=a) <\input| >> anova(reg2) <\output> Analysis of Variance Table \; Response: z \ \ \ \ \ \ \ \ \ \ Df \ Sum Sq Mean Sq \ F value \ \ \ Pr(\F) \ \ \ x \ \ \ \ \ \ \ \ \ 1 255.761 255.761 343.5759 1.706e-14 *** y \ \ \ \ \ \ \ \ \ 1 247.596 247.596 332.6085 2.353e-14 *** x:y \ \ \ \ \ \ \ 1 \ \ 0.063 \ \ 0.063 \ \ 0.0849 \ \ \ 0.7736 \ \ \ Residuals 21 \ 15.633 \ \ 0.744 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> reg3=update(reg2,~.-x:y) <\input| >> anova(reg3) <\output> Analysis of Variance Table \; Response: z \ \ \ \ \ \ \ \ \ \ Df \ Sum Sq Mean Sq F value \ \ \ Pr(\F) \ \ \ x \ \ \ \ \ \ \ \ \ 1 255.761 255.761 \ 358.49 4.174e-15 *** y \ \ \ \ \ \ \ \ \ 1 247.596 247.596 \ 347.04 5.845e-15 *** Residuals 22 \ 15.696 \ \ 0.713 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ <\input| >> \; > both seem significant, so we are left with x and y. <\input| >> reg3 <\output> Call: \ \ \ aov(formula = z ~ x + y, data = a) \; Terms: \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x \ \ \ \ \ \ \ \ y Residuals Sum of Squares \ 255.76053 247.59632 \ 15.69579 Deg. of Freedom \ \ \ \ \ \ \ \ 1 \ \ \ \ \ \ \ \ 1 \ \ \ \ \ \ \ 22 \; Residual standard error: 0.8446567\ Estimated effects may be unbalanced <\input| >> reg3=lm(z~x+y,data=a) <\input| >> reg3 <\output> \; Call: lm(formula = z ~ x + y, data = a) \; Coefficients: (Intercept) \ \ \ \ \ \ \ \ \ \ \ x \ \ \ \ \ \ \ \ \ \ \ y \ \ \ \ \ -0.1195 \ \ \ \ \ -0.9789 \ \ \ \ \ \ 4.9759 \ <\input| >> \; > Now the slope as a function of x is negative! <\input| >> predict(reg3) <\output> \ \ \ \ \ \ \ \ \ \ 1 \ \ \ \ \ \ \ \ \ \ 2 \ \ \ \ \ \ \ \ \ \ 3 \ \ \ \ \ \ \ \ \ \ 4 \ \ \ \ \ \ \ \ \ \ 5 \ \ \ \ \ \ \ \ \ \ 6\ \ 3.87747765 \ 2.89857105 \ 1.91966445 \ 0.94075786 -0.03814874 \ 6.89556959\ \ \ \ \ \ \ \ \ \ \ 7 \ \ \ \ \ \ \ \ \ \ 8 \ \ \ \ \ \ \ \ \ \ 9 \ \ \ \ \ \ \ \ \ 10 \ \ \ \ \ \ \ \ \ 11 \ \ \ \ \ \ \ \ \ 12\ \ 5.91666299 \ 4.93775639 \ 3.95884979 \ 2.97994319 \ 9.91366152 \ 8.93475492\ \ \ \ \ \ \ \ \ \ 13 \ \ \ \ \ \ \ \ \ 14 \ \ \ \ \ \ \ \ \ 15 \ \ \ \ \ \ \ \ \ 16 \ \ \ \ \ \ \ \ \ 17 \ \ \ \ \ \ \ \ \ 18\ \ 7.95584832 \ 6.97694172 \ 5.99803512 12.93175346 11.95284686 10.97394026\ \ \ \ \ \ \ \ \ \ 19 \ \ \ \ \ \ \ \ \ 20 \ \ \ \ \ \ \ \ \ 21 \ \ \ \ \ \ \ \ \ 22 \ \ \ \ \ \ \ \ \ 23 \ \ \ \ \ \ \ \ \ 24\ \ 9.99503366 \ 9.01612706 15.94984539 14.97093879 13.99203219 13.01312559\ \ \ \ \ \ \ \ \ \ 25\ 12.03421899\ <\input| >> plot(z~x,data=a) <\input| >> points(predict(reg3)~x,data=a,col=2);v() <\output> |ps>||||||> <\input| >> \; > \; <\initial> <\collection>