#Source: MapleTrees.R # M <- read.csv("http://www.stats.uwo.ca/faculty/aim/2017/3859/data/MapleTrees.csv") M #View(M) #spreadsheet data viewer str(M) #clearly tree is not numeric! It is a catagorical variable. #we need to make it into a factor variable. TREE <- factor(M$tree, labels=c("Tree 1", "Tree 2", "Tree 3")) M$tree <- TREE head(M) #Oneway ANOVA for differences between trees #The R function to do a one-way anova is aov(). #Note that the function anova() provides an ANOVA summary of an lm-obj. #So aov() and anova() do slightly different things ansOnewayAOV <- aov(velocity~tree, data=M) summary(ansOnewayAOV) # #The lm() function can do anova using dummy variables, ansNull <- lm(velocity~tree, data=M) #eqivalent to one-way ANOVA summary(ansNull) #Note p-value for F-test is same model.matrix(ansNull) #This is the X-matrix aka "design matrix" in regression # #Full Model. Fit separate regression for each tree ansFull <- lm(velocity~tree*loading, data=M) summary(ansFull) model.matrix(ansFull) # #Is the Full Model better than the Null model anova(ansNull, ansFull) #reduced model rejected at < 0.1% # #H0: different intercept but common slope ansCommonSlope <- lm(velocity ~ TREE+loading, data=M) summary(ansCommonSlope) anova(ansCommonSlope) anova(ansCommonSlope, ansFull) #p-value: 5.01%. Accept the common slope model model.matrix(ansCommonSlope) # #diagnostic checks layout(matrix(1:4, ncol=2)) plot(ansCommonSlope, span=1) layout(1) # #residual dependency plot. library(bestglm) res <- resid(ansCommonSlope) X <- model.matrix(ansCommonSlope)[,-1] Xy <- as.data.frame.matrix(cbind(X, res=res)) dgrid(Xy, span=1) # #install.packages("tseries") #for Jarque-Bera test tseries::jarque.bera.test(resid(ansCommonSlope)) #test ok #see CRAN package: nortest, for 5 more omnibus normality tests