<\body> Let us look at the following curious thing: <\input| >> x=rnorm(100) <\input| >> y=rnorm(100,sd=0.1)+x <\input| >> k=ks.test(x,y) <\input| >> k <\output> \; \ \ \ \ \ \ \ \ Two-sample Kolmogorov-Smirnov test \; data: \ x and y\ D = 0.05, p-value = 0.9996 alternative hypothesis: two.sided\ <\input| >> print(k) <\output> \; \ \ \ \ \ \ \ \ Two-sample Kolmogorov-Smirnov test \; data: \ x and y\ D = 0.38, p-value = 1.071e-06 alternative hypothesis: two.sided\ <\input| >> \; > If everything in R is vectors or lists, why is the result of the KS test printed like this? Here is another example: <\input| >> s=summary(x) <\input| >> s <\output> \ \ \ \ Min. \ 1st Qu. \ \ Median \ \ \ \ Mean \ 3rd Qu. \ \ \ \ Max.\ -2.22700 -0.57190 -0.08799 \ 0.01895 \ 0.72390 \ 2.04100\ <\input| >> \; > The magic is the following: Even though things are always vectors or lists, they can have a ''class''. To see the class of an object, we can use the function: <\input| >> class(s) <\output> [1] "table" <\input| >> class(k) <\output> [1] "htest" <\input| >> \; > And now the trick: When we say , the actual function that is called is for k, and for s. <\input| >> print.htest(k) <\output> \; \ \ \ \ \ \ \ \ Two-sample Kolmogorov-Smirnov test \; data: \ x and y\ D = 0.38, p-value = 1.071e-06 alternative hypothesis: two.sided\ <\input| >> print.table(s) <\output> \ \ \ \ Min. \ 1st Qu. \ \ Median \ \ \ \ Mean \ 3rd Qu. \ \ \ \ Max.\ -2.22700 -0.57190 -0.08799 \ 0.01895 \ 0.72390 \ 2.04100\ <\input| >> print.htest(s) <\output> \; Error in strsplit(x, as.character(split), as.logical(extended)) :\ \ \ \ \ \ \ \ \ non-character argument in strsplit() <\input| >> print <\output> function (x, ...)\ UseMethod("print") \environment: namespace:base\ <\input| >> \; > All this is not very important, yet. It only matters if we want to understand what function is called. <\input| >> help(print.htest) <\output> help() for print.htest \ is shown in browser netscape ... Use \ \ \ \ \ help( print.htest , htmlhelp=FALSE) or \ \ \ \ \ \ options(htmlhelp = FALSE) to revert. Warning message:\ Using non-linked HTML file: style sheet and hyperlinks may be incorrect in: help(print.htest)\ <\input| >> \; > Another function that depends on its input is . <\input| >> plot <\output> awk: cmd. line:1: warning: escape sequence `\\.' treated as plain `.' \; function (x, y, ...)\ { \ \ \ \ if (is.null(attr(x, "class")) && is.function(x)) { \ \ \ \ \ \ \ \ nms \- names(list(...)) \ \ \ \ \ \ \ \ if (missing(y))\ \ \ \ \ \ \ \ \ \ \ \ \ y \- { \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if (!"from" %in% nms)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ else if (!"to" %in% nms)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ else if (!"xlim" %in% nms)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ NULL \ \ \ \ \ \ \ \ \ \ \ \ } \ \ \ \ \ \ \ \ if ("ylab" %in% nms)\ \ \ \ \ \ \ \ \ \ \ \ \ plot.function(x, y, ...) \ \ \ \ \ \ \ \ else plot.function(x, y, ylab = paste(deparse(substitute(x)),\ \ \ \ \ \ \ \ \ \ \ \ \ "(x)"), ...) \ \ \ \ } \ \ \ \ else UseMethod("plot") } \environment: namespace:base\ <\input| >> plot(x,y);v() <\output> |ps>||||||> <\input| >> z=as.factor( x\0 ) <\input| >> plot(z,x);v() <\output> |ps>||||||> <\input| >> plot(sin,from=-3,to=3);v() <\output> |ps>||||||> <\input| >> \; > And now we can understand how to bring up help for these functions: <\input| >> class(sin) <\output> [1] "function" <\input| >> help(plot.function) <\output> help() for plot.function \ is shown in browser netscape ... Use \ \ \ \ \ help( plot.function , htmlhelp=FALSE) or \ \ \ \ \ \ options(htmlhelp = FALSE) to revert. Warning message:\ Using non-linked HTML file: style sheet and hyperlinks may be incorrect in: help(plot.function)\ <\input| >> class(z) <\output> awk: cmd. line:1: warning: escape sequence `\\.' treated as plain `.' \; [1] "factor" <\input| >> help(plot.factor) <\output> help() for plot.factor \ is shown in browser netscape ... Use \ \ \ \ \ help( plot.factor , htmlhelp=FALSE) or \ \ \ \ \ \ options(htmlhelp = FALSE) to revert. Warning message:\ Using non-linked HTML file: style sheet and hyperlinks may be incorrect in: help(plot.factor)\ \; > > Last time we learned about qqplot and qqnorm. Let us see how to interpret them: <\input| >> x=rnorm(1000) <\input| >> qqnorm(x);v() <\output> |ps>||||||> <\input| >> abline(0,1,col="red");v() <\output> |ps>||||||> Warning message:\ parameter "color" couldn't be set in high-level plot() function\ <\input| >> \; > abline adds a line to the plot. It has two parameters: slope and intercept. The parameter to plots specifies color. <\input| >> qqnorm(x);abline(0,0,col="red");abline(1,0,col="green");abline(0,1,col="blue") abline(0,2,col="orange");v() <\input| >> \; <\output> |ps>||||||> <\input| >> \; <\input| >> \; > Let us add points on the right to the normal: <\input| >> x=c( rnorm(1000), rnorm(100,mean=-5) ) <\input| >> hist(x,br=30);v() <\output> |ps>||||||> <\input| >> qqnorm(x);abline(0,1,col=2);v() <\output> |ps>||||||> <\input| >> \; > If we have more points on the left than we should have, then quantiles occur earlier, and thus the points are low than they should be. <\input| >> x=rnorm(1000) <\input| >> x=x[ x\-1 ] <\input| >> qqnorm(x); abline(0,1,col=2);v() <\output> |ps>||||||> <\input| >> \; > If we don't have enough points on the left, quantiles happen later, and the points move up. Now let us add points on the right: <\input| >> x=c( rnorm(1000), rnorm(100,mean=3) ) <\input| >> qqnorm(x);abline(0,1,col=2);v() <\output> |ps>||||||> <\input| >> \; > If we have too many points on the right, then the last quantiles occur too late, and the points move up. <\input| >> x=-10:10 ; y= x ^ 2 <\input| >> plot(x,y) ; v() <\output> |ps>||||||> <\input| >> plot(x,y,type="l");v() <\output> |ps>||||||> <\input| >> \; > We can also add to an existing plot <\input| >> plot(x,y);v() <\output> |ps>||||||> <\input| >> points(x,x,col=2);v() <\output> |ps>||||||> <\input| >> abline(40,0,col="green");v() <\output> |ps>||||||> <\input| >> lines(x,y);v() <\output> |ps>||||||> <\input| >> title("A nice plot");v() <\output> |ps>||||||> <\input| >> plot(x,y,pch="M", \ \ \ \ \ col=3, xlab="this is x", ylab= "this is y", \ \ \ \ \ main="A very nice plot" );v() <\input|>> \; <\output> |ps>||||||> \; > The most useful pch is pch=''.'' <\input| >> plot(x,y,pch=".");v() <\output> |ps>||||||> <\input| >> \; > We have a couple of options for saving a plot to a file: One is saving as postscript. For others, look at . First we specify the device, and the file: <\input| >> postscript(file="graph1.eps",width=5,height=5,horizontal=F,onefile=F) <\input| >> \; > Then we do all the plot commands: <\input| >> plot(sin,-3,3) <\input| >> abline(-1,0,col=2); abline(1,0,col=3) <\input| >> \; > Finally we close the device; <\input| >> dev.off() <\output> X11\ \ \ 2\ <\input| >> \; > <\input| >> two.by.two=matrix( 1:4, 2, 2) <\input| >> two.by.two <\output> \ \ \ \ \ [,1] [,2] [1,] \ \ \ 1 \ \ \ 3 [2,] \ \ \ 2 \ \ \ 4 <\input| >> layout(two.by.two) <\input| >> layout.show(4);v() <\output> |ps>||||||> <\input| >> plot(sin,-3,3) <\input| >> plot(cos,-3,3,col=2) <\input| >> plot(tan,-1,1) <\input| >> plot( function(x) x^2, -2, 2) <\input| >> v() <\output> |ps>||||||> <\input| >> layout(1) <\input| >> plot(sin,-3,3);v() <\output> |ps>||||||> <\input| >> l=matrix( c(1,1,2,3),2,2) <\input| >> l <\output> \ \ \ \ \ [,1] [,2] [1,] \ \ \ 1 \ \ \ 2 [2,] \ \ \ 1 \ \ \ 3 <\input| >> layout(l); layout.show(3); v() <\output> |ps>||||||> <\input| >> plot(sin,-3,3); plot(cos,-6,6,col=2); plot(tan,-1.5,1.5,col=3) <\input| >> v() <\output> |ps>||||||> <\input| >> \; > R reserves space for the title, the axes, the labels, and so on. The space is proportional to the character size, and that size does not change when we split the plot. <\input| >> par("mar") <\output> [1] 5.1 4.1 4.1 2.1 <\input| >> \; > above we see how much space is kept for the margin, in lines. top, bottom, left, right. So, if we want to keep figures looking the same, we should scale the text: <\input| >> par(cex=0.5) <\input| >> plot(sin,-3,3); plot(cos,-6,6); plot(tan,-1,1); v() <\output> |ps>||||||> <\input| >> par("cex") <\output> [1] 0.5 <\input| >> layout(1) <\input| >> plot(sin,-3,3); v() <\output> |ps>||||||> <\input| >> par("cex") <\output> [1] 1 <\input| >> \; > The function changes the character size! <\input| >> a=read.table("data/regression.txt",sep="\\t",head=T) <\input| >> a <\output> \ \ growth tannin 1 \ \ \ \ 12 \ \ \ \ \ 0 2 \ \ \ \ 10 \ \ \ \ \ 1 3 \ \ \ \ \ 8 \ \ \ \ \ 2 4 \ \ \ \ 11 \ \ \ \ \ 3 5 \ \ \ \ \ 6 \ \ \ \ \ 4 6 \ \ \ \ \ 7 \ \ \ \ \ 5 7 \ \ \ \ \ 2 \ \ \ \ \ 6 8 \ \ \ \ \ 3 \ \ \ \ \ 7 9 \ \ \ \ \ 3 \ \ \ \ \ 8 <\input| >> plot(a$tannin,a$growth,pch=16);v() <\output> |ps>||||||> <\input| >> \; > The growth of the caterpilar seems to depend on the tannin content. we can try to fit a linear regression. <\input| >> regression=lm( growth ~ tannin, data = a ) <\input| >> regression <\output> \; Call: lm(formula = growth ~ tannin, data = a) \; Coefficients: (Intercept) \ \ \ \ \ \ tannin \ \ \ \ \ \ 11.756 \ \ \ \ \ \ -1.217 \ > this means to do a linear model of growth, so that it depends on tannin. <\input| >> plot( growth ~ tannin, data=a);v() <\output> |ps>||||||> <\input| >> abline( regression );v() <\output> |ps>||||||> <\input| >> \; > R found a good line. The method tries to minimize the square distance to the line. Here is another line: <\input| >> reg1=lm( growth ~ 1, data=a) <\input| >> reg1 <\output> \; Call: lm(formula = growth ~ 1, data = a) \; Coefficients: (Intercept) \ \ \ \ \ \ \ 6.889 \ > In this case, R tries to fit a line that doesn't depend on anything, i.e. it is constant. The best constant to fit the growth is 6.889 <\input| >> abline(6.889,0,col=2); v() <\output> |ps>||||||> <\input| >> mean(a$growth) <\output> [1] 6.888889 <\input| >> \; > Such a line will be a good fit if we have normal errors. For example: <\input| >> x = 1:100 <\input| >> y= 3*x + 15 + rnorm(100,sd=10) <\input| >> plot(x,y); v() <\output> |ps>||||||> <\input| >> reg=lm( y ~ x ) <\input| >> reg <\output> \; Call: lm(formula = y ~ x) \; Coefficients: (Intercept) \ \ \ \ \ \ \ \ \ \ \ x \ \ \ \ \ \ 16.567 \ \ \ \ \ \ \ 2.977 \ <\input| >> \; > We can see that the right values, 3 and 15 were almost found. <\input| >> summary(reg) <\output> \; Call: lm(formula = y ~ x) \; Residuals: \ \ \ \ \ \ Min \ \ \ \ \ \ \ 1Q \ \ \ Median \ \ \ \ \ \ \ 3Q \ \ \ \ \ \ Max\ -31.14652 \ -5.59417 \ -0.03833 \ \ 6.73383 \ 31.15194\ \; Coefficients: \ \ \ \ \ \ \ \ \ \ \ \ Estimate Std. Error t value Pr(\\|t\|) \ \ \ (Intercept) 16.56712 \ \ \ 2.12266 \ \ 7.805 \ 6.6e-12 *** x \ \ \ \ \ \ \ \ \ \ \ 2.97744 \ \ \ 0.03649 \ 81.592 \ \ 2e-16 *** --- Signif. codes: \ 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1\ \; Residual standard error: 10.53 on 98 degrees of freedom Multiple R-Squared: 0.9855, \ \ \ \ Adjusted R-squared: 0.9853\ F-statistic: \ 6657 on 1 and 98 DF, \ p-value: \ 2.2e-16\ <\input| >> \; > Summary tells us that the standard error on the intercept is 2.1, and on the slope 0.03 That means that the estimate of the slope is better than that of the intercept. \; <\initial> <\collection> <\references> <\collection> |?>> |?>> |?>> |?>> |?>> |?>> |?>> |?>>