Motivation
Today we explore various type of regression models and their implementation. First let’s read in some data:
bumpus<-read.csv("Data/bumpus.csv",header=T)
bumpus.data<-log(as.matrix(bumpus[,(5:13)])) # matrix of log-linear measurements
sex<-as.factor(bumpus[,2])
surv<-as.factor(bumpus[,4])
SexBySurv<-as.factor(paste(sex,surv))
Y<-as.matrix(bumpus.data[,1])
X1<-bumpus.data[,2]
X2<-bumpus.data[,3]1: Simple Linear Regression
Simple linear regression fits a bivariate model of the form Y~X, where variation in Y is explained by variation in X. Here are two implementations: parametric and permutational.
lm(Y~X1)##
## Call:
## lm(formula = Y ~ X1)
##
## Coefficients:
## (Intercept) X1
## 1.3045 0.6848#more information
model1<-lm(Y~X1)
summary(model1) #provides regression coefficients ##
## Call:
## lm(formula = Y ~ X1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.034759 -0.010373 -0.002106 0.011899 0.048353
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3045 0.3383 3.856 0.000179 ***
## X1 0.6848 0.0615 11.136 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01617 on 134 degrees of freedom
## Multiple R-squared: 0.4806, Adjusted R-squared: 0.4768
## F-statistic: 124 on 1 and 134 DF, p-value: < 2.2e-16anova(model1) #provides model term tests## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 0.032412 0.032412 124.01 < 2.2e-16 ***
## Residuals 134 0.035023 0.000261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(X1,Y,pch=21, bg="black", cex=2)
abline(model1,lwd=2,col="red")
#Lots of components in model, such as:
model1$fitted.values## 1 2 3 4 5 6 7 8
## 5.060439 5.057592 5.071712 5.091003 5.066099 5.082802 5.071712 5.074501
## 9 10 11 12 13 14 15 16
## 5.077280 5.085547 5.093716 5.077280 5.088281 5.074501 5.077280 5.074501
## 17 18 19 20 21 22 23 24
## 5.077280 5.063275 5.091003 5.074501 5.082802 5.091003 5.085547 5.093716
## 25 26 27 28 29 30 31 32
## 5.066099 5.080046 5.085547 5.091003 5.077280 5.066099 5.093716 5.068911
## 33 34 35 36 37 38 39 40
## 5.077280 5.080046 5.068911 5.093716 5.077280 5.077280 5.080046 5.093716
## 41 42 43 44 45 46 47 48
## 5.054733 5.091003 5.088281 5.088281 5.063275 5.093716 5.088281 5.085547
## 49 50 51 52 53 54 55 56
## 5.077280 5.082802 5.068911 5.046083 5.063275 5.057592 5.071712 5.071712
## 57 58 59 60 61 62 63 64
## 5.068911 5.099108 5.093716 5.051861 5.071712 5.085547 5.077280 5.101788
## 65 66 67 68 69 70 71 72
## 5.077280 5.099108 5.077280 5.074501 5.060439 5.096417 5.074501 5.088281
## 73 74 75 76 77 78 79 80
## 5.057592 5.071712 5.057592 5.080046 5.063275 5.080046 5.034376 5.085547
## 81 82 83 84 85 86 87 88
## 5.085547 5.048978 5.048978 5.093716 5.071712 5.096417 5.071712 5.088281
## 89 90 91 92 93 94 95 96
## 5.031418 5.068911 5.054733 5.077280 5.066099 5.085547 5.071712 5.077280
## 97 98 99 100 101 102 103 104
## 5.077280 5.082802 5.066099 5.066099 5.091003 5.048978 5.028447 5.071712
## 105 106 107 108 109 110 111 112
## 5.057592 5.063275 5.057592 5.051861 5.046083 5.082802 5.066099 5.063275
## 113 114 115 116 117 118 119 120
## 5.048978 5.077280 5.051861 5.051861 5.054733 5.071712 5.043175 5.080046
## 121 122 123 124 125 126 127 128
## 5.051861 5.077280 5.048978 5.057592 5.068911 5.071712 5.074501 5.080046
## 129 130 131 132 133 134 135 136
## 5.071712 5.043175 5.048978 5.068911 5.051861 5.046083 5.074501 5.046083model1$residuals## 1 2 3 4 5
## -2.348683e-02 4.835344e-02 3.461725e-03 -1.582965e-02 -2.267384e-02
## 6 7 8 9 10
## -1.397820e-03 -3.475949e-02 1.309484e-02 -2.742358e-02 8.203329e-03
## 11 12 13 14 15
## -1.231117e-02 1.031675e-02 -3.203480e-02 1.924870e-02 -8.375385e-03
## 16 17 18 19 20
## 6.902868e-03 -1.468455e-02 1.189877e-02 -2.840843e-02 1.309484e-02
## 21 22 23 24 25
## -7.628370e-03 -1.582965e-02 -4.142507e-03 -6.119202e-03 1.530541e-02
## 26 27 28 29 30
## 7.549882e-03 -1.037306e-02 1.494201e-02 -8.375385e-03 1.530541e-02
## 31 32 33 34 35
## -3.112050e-02 1.249309e-02 -8.375385e-03 7.549882e-03 3.095515e-02
## 36 37 38 39 40
## 1.827225e-02 -1.468455e-02 -8.375385e-03 -4.872638e-03 -6.119202e-03
## 41 42 43 44 45
## -4.876740e-03 -9.599096e-03 5.469600e-03 2.370719e-02 -1.341903e-02
## 46 47 48 49 50
## 1.222994e-02 1.766487e-02 2.644092e-02 -2.105772e-03 -7.628370e-03
## 51 52 53 54 55
## -6.316242e-03 3.773460e-03 1.189877e-02 5.002998e-03 -1.546628e-02
## 56 57 58 59 60
## 4.027570e-02 -7.072456e-06 6.837795e-03 -1.854172e-02 4.384334e-03
## 61 62 63 64 65
## -9.117057e-03 1.431956e-02 4.124778e-03 1.019988e-02 -2.105772e-03
## 66 67 68 69 70
## 1.888613e-02 -2.105772e-03 6.902868e-03 -3.000151e-02 1.557087e-02
## 71 72 73 74 75
## -2.464549e-02 -6.876236e-03 -1.416692e-02 -2.185608e-02 -7.736028e-03
## 76 77 78 79 80
## 1.370375e-02 1.189877e-02 1.370375e-02 -1.049579e-02 -1.037306e-02
## 81 82 83 84 85
## -1.037306e-02 8.779018e-04 -5.552989e-03 -6.119202e-03 -1.546628e-02
## 86 87 88 89 90
## -2.666715e-03 3.423338e-02 1.158583e-02 -9.802888e-04 2.483893e-02
## 91 92 93 94 95
## 3.286359e-02 -2.105772e-03 2.149738e-02 -1.037306e-02 -2.807888e-03
## 96 97 98 99 100
## -1.468455e-02 -8.375385e-03 -2.020715e-02 -2.267384e-02 -3.503924e-03
## 101 102 103 104 105
## -3.407126e-03 -5.552989e-03 -4.566759e-03 -2.185608e-02 -2.063943e-02
## 106 107 108 109 110
## 5.629160e-03 -2.715411e-02 -8.436354e-03 -1.564463e-02 1.094802e-02
## 111 112 113 114 115
## -2.267384e-02 3.047516e-02 8.779018e-04 1.647061e-02 4.384334e-03
## 116 117 118 119 120
## 1.704273e-02 -1.130763e-02 9.692275e-03 2.504224e-04 1.981997e-02
## 121 122 123 124 125
## 1.073356e-02 1.031675e-02 -1.854018e-02 5.002998e-03 6.262541e-03
## 126 127 128 129 130
## 1.588425e-02 6.902868e-03 1.981997e-02 -1.546628e-02 1.307111e-02
## 131 132 133 134 135
## 8.779018e-04 -6.316242e-03 -2.142355e-02 -2.657431e-03 1.924870e-02
## 136
## 2.282165e-02#Regression evaluated via residual randomization (RRPP)
library(RRPP)
mydat <- rrpp.data.frame(Y = Y, X1 = X1, X2 = X2, sex = sex, surv = surv, SexBySurv = SexBySurv)
model2 <- lm.rrpp(Y~X1, print.progress = FALSE, data = mydat)
anova(model2)##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Ordinary Least Squares
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## X1 1 0.032412 0.032412 0.48063 124.01 5.9368 0.001 **
## Residuals 134 0.035023 0.000261 0.51937
## Total 135 0.067435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ X1, data = mydat, print.progress = FALSE)anova(model1) #Identical to parametric results ## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 0.032412 0.032412 124.01 < 2.2e-16 ***
## Residuals 134 0.035023 0.000261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coef(model2)## [,1]
## (Intercept) 1.3044593
## X1 0.6847984coef(model1)## (Intercept) X1
## 1.3044593 0.6847984 #plot(model1) #Diagnostic plots: NOTE RUN2: Model II Regression
Model II regresion allows for error in both Y and X, and models deviations orthogonal to the regression line (rather than only in the Y-direction). NOTE: Method is shown for illustrative purposes only. There are very good reasons not to use this for biological data (see lecture for discussion).
library(lmodel2)
lmodel2(Y~X1,nperm=999)## RMA was not requested: it will not be computed.##
## Model II regression
##
## Call: lmodel2(formula = Y ~ X1, nperm = 999)
##
## n = 136 r = 0.6932777 r-square = 0.480634
## Parametric P-values: 2-tailed = 8.511609e-21 1-tailed = 4.255804e-21
## Angle between the two OLS regression lines = 20.53316 degrees
##
## Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
## A permutation test of r is equivalent to a permutation test of the OLS slope
## P-perm for SMA = NA because the SMA slope cannot be tested
##
## Regression results
## Method Intercept Slope Angle (degrees) P-perm (1-tailed)
## 1 OLS 1.3044593 0.6847984 34.40328 0.001
## 2 MA -0.3329159 0.9824064 44.49152 0.001
## 3 SMA -0.3624212 0.9877692 44.64746 NA
##
## Confidence intervals
## Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
## 1 OLS 0.6352913 1.9736273 0.5631720 0.8064248
## 2 MA -1.3892010 0.5532596 0.8213358 1.1743959
## 3 SMA -1.0726264 0.2656984 0.8736027 1.1168556
##
## Eigenvalues: 0.0008563908 0.0001550938
##
## H statistic used for computing C.I. of MA: 0.007883762 RMA<-lmodel2(Y~X1)## RMA was not requested: it will not be computed.## No permutation test will be performedplot(RMA, pch=21,cex=2, bg="black")
abline(model1,lwd=2,col="blue")
3: Multiple Regression
Multiple regression is the case where variation in Y is explained by several independent variables: Y~X1+X2, etc
summary(lm(Y~X1+X2))##
## Call:
## lm(formula = Y ~ X1 + X2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030040 -0.011019 0.000417 0.010957 0.040659
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.82521 0.34843 5.238 6.18e-07 ***
## X1 0.52516 0.07142 7.354 1.77e-11 ***
## X2 0.11042 0.02835 3.895 0.000155 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01537 on 133 degrees of freedom
## Multiple R-squared: 0.5338, Adjusted R-squared: 0.5268
## F-statistic: 76.14 on 2 and 133 DF, p-value: < 2.2e-16anova(lm(Y~X1+X2))## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 0.032412 0.032412 137.118 < 2.2e-16 ***
## X2 1 0.003585 0.003585 15.168 0.0001551 ***
## Residuals 133 0.031438 0.000236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(scatterplot3d)
plot<-scatterplot3d(X1,X2,Y)
plot$plane3d(lm(Y~X1+X2))
One can also perform multiple regression via RRPP. Note that here we have also evaluated the multicollinearity among the independent (X) variables. When this is reasonably large, type II SS may be prefered.
#via RRPP
anova(lm.rrpp(Y~X1+X2,print.progress = FALSE, data=mydat))##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Ordinary Least Squares
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## X1 1 0.032412 0.032412 0.48063 137.118 6.1079 0.001 **
## X2 1 0.003585 0.003585 0.05317 15.168 2.9582 0.001 **
## Residuals 133 0.031438 0.000236 0.46620
## Total 135 0.067435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ X1 + X2, data = mydat, print.progress = FALSE) cor(X1,X2) #hmm, there is multicollinearity in the X-variables. Perhaps use type II SS.## [1] 0.573966anova(lm.rrpp(Y~X1+X2,print.progress = FALSE, data=mydat,SS.type = "II"))##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Ordinary Least Squares
## Sums of Squares and Cross-products: Type II
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## X1 1 0.012782 0.0127819 0.18954 54.074 4.4652 0.001 **
## X2 1 0.003585 0.0035853 0.05317 15.168 2.9582 0.001 **
## Residuals 133 0.031438 0.0002364 0.46620
## Total 135 0.067435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ X1 + X2, SS.type = "II", data = mydat, print.progress = FALSE)A polynomial regression is a special case of multiple regression where X is sequentially raised to higher powers: X, X^2, X^3, etc.
#polynomial regression
fit <- lm(Y~X1) #first degree
fit2 <- lm(Y~poly(X1,2,raw=TRUE))#second degree
fit3 <- lm(Y~poly(X1,3,raw=TRUE))#third degree
fit4 <- lm(Y~poly(X1,4,raw=TRUE))#fourth degree
#evaluate models
anova(fit)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 0.032412 0.032412 124.01 < 2.2e-16 ***
## Residuals 134 0.035023 0.000261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fit,fit2) #In this case, adding polynomial terms NOT an improvement## Analysis of Variance Table
##
## Model 1: Y ~ X1
## Model 2: Y ~ poly(X1, 2, raw = TRUE)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 134 0.035023
## 2 133 0.034977 1 4.6334e-05 0.1762 0.6754anova(fit2,fit3)## Analysis of Variance Table
##
## Model 1: Y ~ poly(X1, 2, raw = TRUE)
## Model 2: Y ~ poly(X1, 3, raw = TRUE)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 133 0.034977
## 2 132 0.034901 1 7.6018e-05 0.2875 0.5927anova(fit3,fit4)## Analysis of Variance Table
##
## Model 1: Y ~ poly(X1, 3, raw = TRUE)
## Model 2: Y ~ poly(X1, 4, raw = TRUE)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 132 0.034901
## 2 132 0.034901 0 0plot(X1,Y)
abline(model1,col="red")
xx <-seq(min(X1),max(X1),length=100)
lines(xx, predict(fit2, data.frame(X1=xx)), col="blue")
lines(xx, predict(fit3, data.frame(X1=xx)), col="green")
lines(xx, predict(fit4, data.frame(X1=xx)), col="purple")## Warning in predict.lm(fit4, data.frame(X1 = xx)): prediction from a rank-
## deficient fit may be misleading
4: ANCOVA: Analysis of Covariance
When one has both categorical and continuous independent X-variables, the model is an ANCOVA. Here one must first test for homogeneity in Y~X slopes across groups. If these are not statistically distinguishable, a common slope model may be used to compare group means while accounting for a common slope.
On the other hand, if the slopes are different, one must perform a comparisons of slopes test to determine how the slopes differ (rather than the group mean test). Examples are below:
anova(lm(Y~X2*SexBySurv))## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X2 1 0.023215 0.0232149 77.3890 8.139e-15 ***
## SexBySurv 3 0.005172 0.0017239 5.7467 0.001014 **
## X2:SexBySurv 3 0.000651 0.0002172 0.7239 0.539496
## Residuals 128 0.038397 0.0003000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Implies slopes for M and F statistically equivalent, so drop and re-run
anova(lm(Y~X2+SexBySurv))## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X2 1 0.023215 0.0232149 77.8814 6.015e-15 ***
## SexBySurv 3 0.005172 0.0017239 5.7833 0.0009591 ***
## Residuals 131 0.039048 0.0002981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model.ancova<-lm(Y~X2+SexBySurv)
colo<-rep("blue",nrow(bumpus.data)); colo[which(SexBySurv == 'f TRUE')]<-"red";
colo[which(SexBySurv == 'f FALSE')]<-"pink"; colo[which(SexBySurv == 'm FALSE')]<-"lightblue";
plot(X2,Y,pch=21,cex=2,bg=colo)
abline(model.ancova$coefficients[1],model.ancova$coefficients[2])
abline((model.ancova$coefficients[1]+model.ancova$coefficients[3]),model.ancova$coefficients[2])
abline((model.ancova$coefficients[1]+model.ancova$coefficients[4]),model.ancova$coefficients[2])
abline((model.ancova$coefficients[1]+model.ancova$coefficients[5]),model.ancova$coefficients[2])
#note: parameters are intercept, slope, deviations of int. for groups
#via RRPP
model.anc2 <- lm.rrpp(Y~X2*SexBySurv, print.progress = FALSE, data = mydat)
anova(model.anc2)##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Ordinary Least Squares
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## X2 1 0.023215 0.0232149 0.34426 77.3890 5.2846 0.001 **
## SexBySurv 3 0.005172 0.0017239 0.07669 5.7467 3.0048 0.002 **
## X2:SexBySurv 3 0.000651 0.0002172 0.00966 0.7239 -0.0862 0.530
## Residuals 128 0.038397 0.0003000 0.56939
## Total 135 0.067435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ X2 * SexBySurv, data = mydat, print.progress = FALSE)model.anc3 <- lm.rrpp(Y~X2+SexBySurv, print.progress = FALSE, data = mydat)
#NOTE: `plot` works on lm.rrpp objects, but is most useful for #multivariate data
plot(X2,Y,pch=21,cex=2,bg=colo)
abline(model.ancova$coefficients[1],model.ancova$coefficients[2])
abline((model.ancova$coefficients[1]+model.ancova$coefficients[3]),model.ancova$coefficients[2])
abline((model.ancova$coefficients[1]+model.ancova$coefficients[4]),model.ancova$coefficients[2])
abline((model.ancova$coefficients[1]+model.ancova$coefficients[5]),model.ancova$coefficients[2])
With ANCOVA, one can still perform pairwise comparisons among groups. However, which comparisons to perform depends upon the outcome of the ANCOVA. For example, if the test of slopes model (i.e., interaction term) did not reveal differences, then a common slopes model is fit to the data. In this case, one compares differences among groups in terms of their least squares means. On the other hand, if the interaction term implied that one or more groups display differing slopes, then pairwise comparisons of slopes are in order. The former may be accomplished using ‘standard’ approaches. Both tests of means and tests of slopes are facilitated in RRPP (see help files for publications describing the theoretical developments of these methods).
### Pairwise comparisons: test of means
pairwise.t.test(model.ancova$fitted.values, SexBySurv, p.adj = "none") #standard approach##
## Pairwise comparisons using t tests with pooled SD
##
## data: model.ancova$fitted.values and SexBySurv
##
## f FALSE f TRUE m FALSE
## f TRUE 0.028 - -
## m FALSE 4.2e-15 < 2e-16 -
## m TRUE 0.029 1.6e-05 1.0e-12
##
## P value adjustment method: none#Pairwise comparisons of means via RRPP
summary(pairwise(model.anc3, groups = SexBySurv), test.type = "dist", stat.table = FALSE) ##
## Pairwise comparisons
##
## Groups: f FALSE f TRUE m FALSE m TRUE
##
## RRPP: 1000 permutations
##
## LS means:
## Vectors hidden (use show.vectors = TRUE to view)
##
## Pairwise distances between means
## f FALSE f TRUE m FALSE m TRUE
## f FALSE 0.000000000 0.001283828 0.01595601 0.004128217
## f TRUE 0.001283828 0.000000000 0.01723984 0.005412045
## m FALSE 0.015956012 0.017239840 0.00000000 0.011827794
## m TRUE 0.004128217 0.005412045 0.01182779 0.000000000
##
## Pairwise 95% Upper confidence limits between means
## f FALSE f TRUE m FALSE m TRUE
## f FALSE 0.000000000 0.009669733 0.008780365 0.008878357
## f TRUE 0.009669733 0.000000000 0.010341481 0.009216073
## m FALSE 0.008780365 0.010341481 0.000000000 0.008026951
## m TRUE 0.008878357 0.009216073 0.008026951 0.000000000
##
## Pairwise effect sizes (Z) between means
## f FALSE f TRUE m FALSE m TRUE
## f FALSE 0.0000000 -0.9612166 2.772571 0.4670408
## f TRUE -0.9612166 0.0000000 2.578038 0.6999764
## m FALSE 2.7725709 2.5780381 0.000000 2.3954073
## m TRUE 0.4670408 0.6999764 2.395407 0.0000000
##
## Pairwise P-values between means
## f FALSE f TRUE m FALSE m TRUE
## f FALSE 1.000 0.819 0.002 0.337
## f TRUE 0.819 1.000 0.002 0.256
## m FALSE 0.002 0.002 1.000 0.002
## m TRUE 0.337 0.256 0.002 1.000#Pairwise comparisons of SLOPES (NOTE: makes sense only when interaction term significant)
PW <- pairwise(model.anc3, groups = SexBySurv)
summary(PW, test.type = "VC", angle.type = "deg") #comparison of angles between vectors##
## Pairwise comparisons
##
## Groups: f FALSE f TRUE m FALSE m TRUE
##
## RRPP: 1000 permutations
##
## LS means:
## Vectors hidden (use show.vectors = TRUE to view)
##
## Pairwise statistics based on mean vector correlations
## r angle UCL (95%) Z Pr > angle
## f FALSE:f TRUE 1 0 8.537736e-07 -0.4834673 0.5960
## f FALSE:m FALSE 1 0 8.537736e-07 -0.4760313 0.5935
## f FALSE:m TRUE 1 0 8.537736e-07 -0.4740073 0.5930
## f TRUE:m FALSE 1 0 8.537736e-07 -0.4781134 0.5945
## f TRUE:m TRUE 1 0 8.537736e-07 -0.4832523 0.5960
## m FALSE:m TRUE 1 0 8.537736e-07 -0.4830427 0.59555: Permutation Test for Regression
Here is code for performing a permutation test to evaluate linear regression. It is written in ‘loop’ style, so it is more readable (in a few weeks we’ll discuss other implementations). Note that the function ‘procD.lm’ in the package ‘geomorph’ also performs this test.
F.obs<-anova(lm(Y~X1))$F[[1]] #Find Test value and save
permute<-1999
F.rand.vec<-array(NA,(permute+1))
F.rand.vec[permute+1]<-F.obs
for(i in 1:permute){
###Shuffle Data
y.rand<-sample(Y) #Resample vector
F.rand.vec[i]<-anova(lm(y.rand~X1))$F[[1]]
}
F.obs## [1] 124.0069P.Ftest<-rank(F.rand.vec[permute+1])/(permute+1)
P.Ftest## [1] 5e-04####Plot
hist(F.rand.vec,40,freq=T,col="gray")
segments(F.obs, 0, F.obs, 50) ##Plot Observed value