sires=read.table("siredata.txt", header=T,sep="\t",skip=3)
prog=read.table("progdata.txt", header=T,sep="\t",skip=3)
dim(sires)
## [1] 10 12
dim(prog)
## [1] 400 14
print(sires)
## id weight m11 m12 m21 m22 m31 m32 m41 m42 m51 m52
## 1 sire1 334.14 M2 M1 M3 M2 M3 M4 M4 M2 M4 M2
## 2 sire2 364.81 M3 M2 M3 M2 M2 M3 M2 M4 M3 M1
## 3 sire3 383.95 M2 M4 M2 M4 M3 M2 M3 M4 M1 M4
## 4 sire4 349.88 M2 M1 M1 M2 M4 M3 M2 M1 M4 M3
## 5 sire5 357.87 M1 M3 M2 M1 M3 M1 M3 M4 M3 M2
## 6 sire6 364.87 M3 M2 M1 M3 M4 M2 M3 M1 M1 M3
## 7 sire7 361.36 M4 M1 M4 M1 M2 M1 M1 M2 M2 M4
## 8 sire8 357.56 M3 M4 M2 M4 M4 M1 M2 M4 M1 M3
## 9 sire9 333.49 M1 M2 M4 M3 M3 M2 M1 M2 M3 M2
## 10 sire10 360.92 M2 M3 M3 M1 M3 M2 M2 M3 M3 M2
head(prog)
## id sire sex weight m11 m12 m21 m22 m31 m32 m41 m42 m51 m52
## 1 id1 sire1 F 293.61 M5 M6 M2 M5 M3 M6 M2 M5 M4 M4
## 2 id2 sire1 M 335.43 M1 M4 M3 M1 M4 M5 M4 M3 M2 M2
## 3 id3 sire1 M 340.09 M2 M3 M2 M6 M3 M1 M4 M3 M4 M3
## 4 id4 sire1 M 343.08 M2 M3 M2 M1 M4 M6 M2 M3 M4 M5
## 5 id5 sire1 F 287.08 M1 M3 M2 M4 M4 M6 M4 M5 M2 M3
## 6 id6 sire1 F 302.17 M2 M2 M2 M5 M3 M5 M2 M4 M4 M1
summary(sires)
## id weight m11 m12 m21 m22 m31 m32
## sire1 :1 Min. :333.5 M1:2 M1:3 M1:2 M1:3 M2:2 M1:3
## sire10 :1 1st Qu.:351.8 M2:4 M2:3 M2:3 M2:3 M3:5 M2:4
## sire2 :1 Median :359.4 M3:3 M3:2 M3:3 M3:2 M4:3 M3:2
## sire3 :1 Mean :356.9 M4:1 M4:2 M4:2 M4:2 M4:1
## sire4 :1 3rd Qu.:363.9
## sire5 :1 Max. :383.9
## (Other):4
## m41 m42 m51 m52
## M1:2 M1:2 M1:3 M1:1
## M2:4 M2:3 M2:1 M2:4
## M3:3 M3:1 M3:4 M3:3
## M4:1 M4:4 M4:2 M4:2
##
##
##
summary(prog)
## id sire sex weight m11 m12
## id1 : 1 sire1 : 40 F:200 304.06 : 3 - : 1 - : 1
## id10 : 1 sire10 : 40 M:200 299.28 : 2 M1: 98 M1:70
## id100 : 1 sire2 : 40 307.21 : 2 M2:139 M2:68
## id101 : 1 sire3 : 40 311.7 : 2 M3:101 M3:67
## id102 : 1 sire4 : 40 313.62 : 2 M4: 60 M4:62
## id103 : 1 sire5 : 40 313.66 : 2 M5: 1 M5:74
## (Other):394 (Other):160 (Other):387 M6:58
## m21 m22 m31 m32 m41 m42 m51 m52
## - : 1 - : 1 - : 1 - : 1 - : 1 - : 1 - : 1 - : 1
## M1: 97 M1:64 M1: 55 M1:66 M1: 78 M1:57 M1: 83 M1:58
## M2:130 M2:63 M2:115 M2:69 M2:142 M2:68 M2: 94 M2:64
## M3: 93 M3:79 M3:134 M3:72 M3: 82 M3:78 M3:135 M3:64
## M4: 79 M4:61 M4: 95 M4:61 M4: 97 M4:75 M4: 87 M4:74
## M5:56 M5:69 M5:68 M5:76
## M6:76 M6:62 M6:53 M6:63
prog[which(prog$weight=="-"),]
## id sire sex weight m11 m12 m21 m22 m31 m32 m41 m42 m51 m52
## 90 id90 sire3 M - M4 M1 M2 M5 M2 M5 M4 M2 M1 M1
prog=prog[-which(prog$weight=="-"),]
prog$weight=as.numeric(as.character(prog$weight))
summary(prog$weight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 282.8 307.0 330.9 329.9 351.8 480.5
summary(sires$weight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 333.5 351.8 359.4 356.9 363.9 383.9
plot(prog$weight,main="XY plot of weight by animal", xlab="animal",ylab="weight",col="green")
prog=prog[-which(prog$weight>400),]
summary(prog$weight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 282.8 307.0 330.1 329.5 351.6 386.3
index=grep("m",names(prog))
missing=numeric()
for (i in 1:length(index)) missing=c(missing,which(prog[,index[i]]=="-"))
print(missing)
## [1] 69 69 69 69 69 69 69 69 69 69
missingU=unique(missing)
print(missingU)
## [1] 69
print(prog[missingU,])
## id sire sex weight m11 m12 m21 m22 m31 m32 m41 m42 m51 m52
## 70 id70 sire2 M 372.45 - - - - - - - - - -
prog=prog[-missingU,]
summary(prog$m11)
## - M1 M2 M3 M4 M5
## 0 97 139 101 59 1
for (i in 1:length(index))
prog[,index[i]]=factor(prog[,index[i]])
summary(prog$m11)
## M1 M2 M3 M4 M5
## 97 139 101 59 1
boxplot(prog$weight~prog$sire, col=1:length(levels(prog$sire)),
main="Boxplot of weights by sire")
boxplot(prog$weight~prog$sex, col=1:length(levels(prog$sex)),
main="Boxplot of weights by sex")
plot(density(prog$weight),col="blue", main="Density plot of weights")
plot(prog$weight,col=prog$sex, pch=as.numeric (prog$sex),main="XY plot of weight by animal", xlab="animal",ylab="weight")
legend("topleft",levels(prog$sex),col=1:2,pch=1:2)
indexms=grep("m",names(sires))
indexms=matrix(indexms,length(indexms)/2,2,byrow=T)
print(indexms)
## [,1] [,2]
## [1,] 3 4
## [2,] 5 6
## [3,] 7 8
## [4,] 9 10
## [5,] 11 12
indexm=grep("m",names(prog))
indexm=matrix(indexm,length(indexm)/2,2,byrow=T)
print(indexm)
## [,1] [,2]
## [1,] 5 6
## [2,] 7 8
## [3,] 9 10
## [4,] 11 12
## [5,] 13 14
alleles=summary(factor(c(as.character(prog$m11), as.character(prog$m12))))
print(alleles)
## M1 M2 M3 M4 M5 M6
## 166 207 167 121 75 58
alleles=alleles/sum(alleles)
print(alleles)
## M1 M2 M3 M4 M5 M6
## 0.20906801 0.26070529 0.21032746 0.15239295 0.09445844 0.07304786
hold=data.frame(m11=as.character(prog$m11), m12=as.character(prog$m12))
hold[,1]=as.character(hold[,1])
hold[,2]=as.character(hold[,2])
sorted=character()
for (i in 1:length(hold[,1]))
sorted=rbind(sorted,sort(as.character(hold[i,])))
genotypes=paste(as.character(sorted[,1]), as.character(sorted[,2]),sep="_")
genotypes=summary(factor(genotypes))
print(genotypes)
## M1_M1 M1_M2 M1_M3 M1_M4 M1_M5 M1_M6 M2_M2 M2_M3 M2_M4 M2_M5 M2_M6 M3_M3
## 14 42 34 22 26 14 28 35 38 17 19 21
## M3_M4 M3_M5 M3_M6 M4_M4 M4_M5 M4_M6 M5_M6
## 24 19 13 7 12 11 1
genotypes=genotypes/sum(genotypes)
print(genotypes)
## M1_M1 M1_M2 M1_M3 M1_M4 M1_M5 M1_M6
## 0.035264484 0.105793451 0.085642317 0.055415617 0.065491184 0.035264484
## M2_M2 M2_M3 M2_M4 M2_M5 M2_M6 M3_M3
## 0.070528967 0.088161209 0.095717884 0.042821159 0.047858942 0.052896725
## M3_M4 M3_M5 M3_M6 M4_M4 M4_M5 M4_M6
## 0.060453401 0.047858942 0.032745592 0.017632242 0.030226700 0.027707809
## M5_M6
## 0.002518892
barplot(sort(alleles,decreasing=T),col=1:11, main="Barplot of allelic frequencies")
barplot(sort(genotypes,decreasing=T),col=1:11, main="Barplot of genotypic frequencies")
allgeno=NULL
for(i in 1:length(indexm[,1]))
{
hold=data.frame(prog[indexm[i,]])
hold[,1]=as.character(hold[,1])
hold[,2]=as.character(hold[,2])
sorted=character()
for (i in 1:length(hold[,1]))
sorted=rbind(sorted,sort(as.character(hold[i,])))
genotypes=paste(as.character(sorted[,1]), as.character(sorted[,2]),sep="_")
allgeno=cbind(allgeno,genotypes)
}
colnames(allgeno)=c("M1","M2","M3","M4","M5")
markers=data.frame(prog[,1:4],allgeno)
head(markers)
## id sire sex weight M1 M2 M3 M4 M5
## 1 id1 sire1 F 293.61 M5_M6 M2_M5 M3_M6 M2_M5 M4_M4
## 2 id2 sire1 M 335.43 M1_M4 M1_M3 M4_M5 M3_M4 M2_M2
## 3 id3 sire1 M 340.09 M2_M3 M2_M6 M1_M3 M3_M4 M3_M4
## 4 id4 sire1 M 343.08 M2_M3 M1_M2 M4_M6 M2_M3 M4_M5
## 5 id5 sire1 F 287.08 M1_M3 M2_M4 M4_M6 M4_M5 M2_M3
## 6 id6 sire1 F 302.17 M2_M2 M2_M5 M3_M5 M2_M4 M1_M4
lm(weight~M1,data=markers)
##
## Call:
## lm(formula = weight ~ M1, data = markers)
##
## Coefficients:
## (Intercept) M1M1_M2 M1M1_M3 M1M1_M4 M1M1_M5
## 3.253e+02 1.771e+00 2.256e+00 8.919e+00 -1.372e+00
## M1M1_M6 M1M2_M2 M1M2_M3 M1M2_M4 M1M2_M5
## 7.143e-04 -1.049e+00 9.872e+00 -3.144e-01 -2.640e+00
## M1M2_M6 M1M3_M3 M1M3_M4 M1M3_M5 M1M3_M6
## 1.145e+01 4.787e+00 7.156e+00 1.295e+01 8.932e+00
## M1M4_M4 M1M4_M5 M1M4_M6 M1M5_M6
## 1.955e+01 1.065e+01 6.623e-02 -3.164e+01
results=lm(weight~M1,data=markers)
class(results)
## [1] "lm"
summary(results)
##
## Call:
## lm(formula = weight ~ M1, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.185 -22.034 1.259 21.156 52.109
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.253e+02 6.763e+00 48.095 <2e-16 ***
## M1M1_M2 1.771e+00 7.809e+00 0.227 0.8207
## M1M1_M3 2.256e+00 8.035e+00 0.281 0.7790
## M1M1_M4 8.919e+00 8.651e+00 1.031 0.3032
## M1M1_M5 -1.372e+00 8.388e+00 -0.164 0.8702
## M1M1_M6 7.143e-04 9.564e+00 0.000 0.9999
## M1M2_M2 -1.049e+00 8.283e+00 -0.127 0.8993
## M1M2_M3 9.872e+00 8.002e+00 1.234 0.2181
## M1M2_M4 -3.144e-01 7.911e+00 -0.040 0.9683
## M1M2_M5 -2.640e+00 9.132e+00 -0.289 0.7727
## M1M2_M6 1.145e+01 8.913e+00 1.284 0.1998
## M1M3_M3 4.787e+00 8.731e+00 0.548 0.5838
## M1M3_M4 7.156e+00 8.510e+00 0.841 0.4009
## M1M3_M5 1.295e+01 8.913e+00 1.453 0.1472
## M1M3_M6 8.932e+00 9.746e+00 0.916 0.3600
## M1M4_M4 1.955e+01 1.171e+01 1.669 0.0959 .
## M1M4_M5 1.065e+01 9.954e+00 1.069 0.2856
## M1M4_M6 6.623e-02 1.020e+01 0.006 0.9948
## M1M5_M6 -3.164e+01 2.619e+01 -1.208 0.2278
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.3 on 378 degrees of freedom
## Multiple R-squared: 0.04804, Adjusted R-squared: 0.00271
## F-statistic: 1.06 on 18 and 378 DF, p-value: 0.3918
plot(predict(results),residuals(results), main="XY plot of residuals X fitted values",
xlab="fitted values (weight)",
ylab="residuals",col="blue")
qqnorm(predict(results),col="blue")
qqnorm(residuals(results),col="blue",main="")
fligner.test(weight~M1,data=markers)
##
## Fligner-Killeen test of homogeneity of variances
##
## data: weight by M1
## Fligner-Killeen:med chi-squared = 9.398, df = 18, p-value = 0.9498
anova(results)
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## M1 18 12214 678.56 1.0598 0.3918
## Residuals 378 242028 640.29
summary(lm(weight~sex+M1,data=markers))
##
## Call:
## lm(formula = weight ~ sex + M1, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.1174 -7.2282 -0.7501 7.1186 31.3916
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.992e+02 2.900e+00 103.168 < 2e-16 ***
## sexM 4.558e+01 1.081e+00 42.143 < 2e-16 ***
## M1M1_M2 6.111e+00 3.274e+00 1.867 0.062693 .
## M1M1_M3 5.512e+00 3.368e+00 1.637 0.102537
## M1M1_M4 1.010e+01 3.625e+00 2.787 0.005590 **
## M1M1_M5 1.884e+00 3.516e+00 0.536 0.592370
## M1M1_M6 7.143e-04 4.007e+00 0.000 0.999858
## M1M2_M2 7.090e+00 3.476e+00 2.040 0.042079 *
## M1M2_M3 5.965e+00 3.354e+00 1.779 0.076118 .
## M1M2_M4 7.739e+00 3.320e+00 2.331 0.020292 *
## M1M2_M5 4.637e+00 3.830e+00 1.211 0.226794
## M1M2_M6 1.350e+01 3.735e+00 3.615 0.000341 ***
## M1M3_M3 9.128e+00 3.660e+00 2.494 0.013051 *
## M1M3_M4 1.421e+01 3.569e+00 3.981 8.24e-05 ***
## M1M3_M5 1.260e+01 3.734e+00 3.375 0.000815 ***
## M1M3_M6 1.043e+01 4.084e+00 2.555 0.011010 *
## M1M4_M4 1.304e+01 4.910e+00 2.655 0.008258 **
## M1M4_M5 1.010e+01 4.171e+00 2.422 0.015898 *
## M1M4_M6 9.537e+00 4.278e+00 2.229 0.026375 *
## M1M5_M6 -5.599e+00 1.099e+01 -0.509 0.610810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.6 on 377 degrees of freedom
## Multiple R-squared: 0.8333, Adjusted R-squared: 0.8249
## F-statistic: 99.2 on 19 and 377 DF, p-value: < 2.2e-16
summary(lm(weight~sex+sire+M1,data=markers))
##
## Call:
## lm(formula = weight ~ sex + sire + M1, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.5650 -5.6855 -0.0582 5.2260 21.5778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 293.8628 2.6104 112.574 < 2e-16 ***
## sexM 45.5573 0.8566 53.186 < 2e-16 ***
## siresire10 14.5459 1.9541 7.444 7.00e-13 ***
## siresire2 18.2750 2.0120 9.083 < 2e-16 ***
## siresire3 27.9691 2.0110 13.908 < 2e-16 ***
## siresire4 8.5190 1.9165 4.445 1.16e-05 ***
## siresire5 14.2383 1.9915 7.150 4.72e-12 ***
## siresire6 15.8241 2.0043 7.895 3.37e-14 ***
## siresire7 14.2136 2.0673 6.875 2.66e-11 ***
## siresire8 14.1905 2.1265 6.673 9.22e-11 ***
## siresire9 3.6338 1.9234 1.889 0.0596 .
## M1M1_M2 1.0713 2.6239 0.408 0.6833
## M1M1_M3 -0.7715 2.7549 -0.280 0.7796
## M1M1_M4 0.9751 3.0089 0.324 0.7461
## M1M1_M5 -1.6311 2.8036 -0.582 0.5611
## M1M1_M6 -1.6844 3.1542 -0.534 0.5937
## M1M2_M2 -0.4491 2.8385 -0.158 0.8744
## M1M2_M3 -2.0030 2.7399 -0.731 0.4652
## M1M2_M4 -1.4636 2.7087 -0.540 0.5893
## M1M2_M5 -1.9519 3.0634 -0.637 0.5244
## M1M2_M6 2.4693 3.0634 0.806 0.4207
## M1M3_M3 -0.8174 3.0775 -0.266 0.7907
## M1M3_M4 2.1514 2.9849 0.721 0.4715
## M1M3_M5 2.2666 3.1662 0.716 0.4745
## M1M3_M6 0.4037 3.3847 0.119 0.9051
## M1M4_M4 -3.6682 4.0730 -0.901 0.3684
## M1M4_M5 1.2526 3.6037 0.348 0.7283
## M1M4_M6 -3.0625 3.6157 -0.847 0.3975
## M1M5_M6 -0.2528 8.7037 -0.029 0.9768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.303 on 368 degrees of freedom
## Multiple R-squared: 0.9002, Adjusted R-squared: 0.8926
## F-statistic: 118.6 on 28 and 368 DF, p-value: < 2.2e-16
model1=lm(weight~sex+sire+M1,data=markers)
model2=lm(weight~sex+M1,data=markers)
model3=lm(weight~M1,data=markers)
anova(model1,model2)
## Analysis of Variance Table
##
## Model 1: weight ~ sex + sire + M1
## Model 2: weight ~ sex + M1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 368 25370
## 2 377 42379 -9 -17009 27.413 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model2,model3)
## Analysis of Variance Table
##
## Model 1: weight ~ sex + M1
## Model 2: weight ~ M1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 377 42379
## 2 378 242028 -1 -199649 1776.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model4=lm(weight~1,data=markers)
anova(model3,model4)
## Analysis of Variance Table
##
## Model 1: weight ~ M1
## Model 2: weight ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 378 242028
## 2 396 254242 -18 -12214 1.0598 0.3918
anova(model2,model4)
## Analysis of Variance Table
##
## Model 1: weight ~ sex + M1
## Model 2: weight ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 377 42379
## 2 396 254242 -19 -211863 99.196 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(weight~sex+M1, data=markers[which(prog$sire=="sire1"),]))
##
## Call:
## lm(formula = weight ~ sex + M1, data = markers[which(prog$sire ==
## "sire1"), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.196 -3.861 0.000 3.622 12.664
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 292.0005 3.9627 73.687 <2e-16 ***
## sexM 44.0293 2.4834 17.730 <2e-16 ***
## M1M1_M2 4.0354 4.1938 0.962 0.344
## M1M1_M3 -0.7563 4.6992 -0.161 0.873
## M1M1_M4 -0.5998 7.2481 -0.083 0.935
## M1M1_M5 2.6276 4.7673 0.551 0.586
## M1M1_M6 4.2302 7.2481 0.584 0.564
## M1M2_M2 1.7728 5.3540 0.331 0.743
## M1M2_M3 5.1777 4.8342 1.071 0.294
## M1M2_M4 1.0178 4.6020 0.221 0.827
## M1M2_M5 -0.1500 5.0916 -0.029 0.977
## M1M5_M6 1.6095 7.3885 0.218 0.829
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.236 on 27 degrees of freedom
## Multiple R-squared: 0.9497, Adjusted R-squared: 0.9292
## F-statistic: 46.35 on 11 and 27 DF, p-value: 1.26e-14
summary(lm(weight~sex+sire+M2,data=markers))
##
## Call:
## lm(formula = weight ~ sex + sire + M2, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.7861 -5.6922 -0.5459 4.9430 23.5192
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 294.02282 2.54891 115.352 < 2e-16 ***
## sexM 45.62156 0.85799 53.172 < 2e-16 ***
## siresire10 16.00473 1.94028 8.249 2.86e-15 ***
## siresire2 19.39024 1.92604 10.067 < 2e-16 ***
## siresire3 27.62483 2.05790 13.424 < 2e-16 ***
## siresire4 9.09786 1.99472 4.561 6.94e-06 ***
## siresire5 14.96327 1.97928 7.560 3.23e-13 ***
## siresire6 17.28381 1.95646 8.834 < 2e-16 ***
## siresire7 15.56163 2.07181 7.511 4.46e-13 ***
## siresire8 14.58313 2.02544 7.200 3.40e-12 ***
## siresire9 3.68223 1.96666 1.872 0.0620 .
## M2M1_M2 -1.35057 2.51831 -0.536 0.5921
## M2M1_M3 -1.90570 2.40153 -0.794 0.4280
## M2M1_M4 -4.73343 2.72489 -1.737 0.0832 .
## M2M1_M5 -5.52770 2.96021 -1.867 0.0626 .
## M2M1_M6 0.07987 2.73991 0.029 0.9768
## M2M2_M2 -1.03442 2.71376 -0.381 0.7033
## M2M2_M3 -1.36150 2.49168 -0.546 0.5851
## M2M2_M4 -0.64962 2.51417 -0.258 0.7963
## M2M2_M5 1.20683 2.81852 0.428 0.6688
## M2M2_M6 -1.03355 2.65616 -0.389 0.6974
## M2M3_M3 -1.88421 2.78631 -0.676 0.4993
## M2M3_M4 -0.71071 2.67845 -0.265 0.7909
## M2M3_M5 -4.11478 3.12493 -1.317 0.1887
## M2M3_M6 2.29599 2.93780 0.782 0.4350
## M2M4_M4 0.17240 2.94429 0.059 0.9533
## M2M4_M5 -0.06417 3.82402 -0.017 0.9866
## M2M4_M6 0.01216 3.13463 0.004 0.9969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.274 on 369 degrees of freedom
## Multiple R-squared: 0.9006, Adjusted R-squared: 0.8934
## F-statistic: 123.9 on 27 and 369 DF, p-value: < 2.2e-16
summary(lm(weight~sex+sire+M3,data=markers))
##
## Call:
## lm(formula = weight ~ sex + sire + M3, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.3468 -5.8351 -0.3536 4.6471 23.8249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 295.7760 3.3993 87.011 < 2e-16 ***
## sexM 45.3631 0.8542 53.107 < 2e-16 ***
## siresire10 14.5490 2.0061 7.252 2.43e-12 ***
## siresire2 19.8474 1.9860 9.993 < 2e-16 ***
## siresire3 27.9468 1.9841 14.085 < 2e-16 ***
## siresire4 8.6667 1.9207 4.512 8.64e-06 ***
## siresire5 15.1702 2.0052 7.565 3.11e-13 ***
## siresire6 16.4811 1.9533 8.438 7.45e-16 ***
## siresire7 15.0670 2.0644 7.299 1.80e-12 ***
## siresire8 14.5860 1.9627 7.432 7.54e-13 ***
## siresire9 3.2690 1.9850 1.647 0.100
## M3M1_M2 -1.7042 3.4473 -0.494 0.621
## M3M1_M3 -2.7984 3.3134 -0.845 0.399
## M3M1_M4 -3.5918 3.5019 -1.026 0.306
## M3M1_M5 -5.2671 3.7465 -1.406 0.161
## M3M1_M6 -2.9972 3.9576 -0.757 0.449
## M3M2_M2 -3.5300 3.7035 -0.953 0.341
## M3M2_M3 -1.7404 3.3256 -0.523 0.601
## M3M2_M4 -3.1059 3.3188 -0.936 0.350
## M3M2_M5 -5.2810 3.5598 -1.484 0.139
## M3M2_M6 2.7602 3.9160 0.705 0.481
## M3M3_M3 -3.6036 3.6003 -1.001 0.318
## M3M3_M4 -3.4700 3.3350 -1.040 0.299
## M3M3_M5 -1.1775 3.6049 -0.327 0.744
## M3M3_M6 -1.9506 3.5155 -0.555 0.579
## M3M4_M4 -2.0059 3.8403 -0.522 0.602
## M3M4_M5 0.4208 3.8969 0.108 0.914
## M3M4_M6 -1.4888 3.7960 -0.392 0.695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.292 on 369 degrees of freedom
## Multiple R-squared: 0.9002, Adjusted R-squared: 0.8929
## F-statistic: 123.3 on 27 and 369 DF, p-value: < 2.2e-16
summary(lm(weight~sex+sire+M4,data=markers))
##
## Call:
## lm(formula = weight ~ sex + sire + M4, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.7807 -5.8044 -0.0684 5.2611 23.2223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 295.6114 3.1535 93.742 < 2e-16 ***
## sexM 45.7380 0.8573 53.354 < 2e-16 ***
## siresire10 15.7838 1.9385 8.142 6.04e-15 ***
## siresire2 19.7070 1.9098 10.319 < 2e-16 ***
## siresire3 28.1573 1.9728 14.273 < 2e-16 ***
## siresire4 7.6642 2.0358 3.765 0.000194 ***
## siresire5 14.7311 1.9564 7.530 3.94e-13 ***
## siresire6 16.2489 2.0596 7.889 3.48e-14 ***
## siresire7 14.5967 1.9746 7.392 9.77e-13 ***
## siresire8 15.1843 1.9231 7.896 3.33e-14 ***
## siresire9 3.1462 1.9987 1.574 0.116322
## M4M1_M2 -4.0779 2.9287 -1.392 0.164631
## M4M1_M3 -3.2467 3.1811 -1.021 0.308106
## M4M1_M4 -2.2062 3.0829 -0.716 0.474682
## M4M1_M5 -0.2962 3.4820 -0.085 0.932263
## M4M1_M6 5.3292 3.8048 1.401 0.162154
## M4M2_M2 -3.2303 3.3115 -0.975 0.329967
## M4M2_M3 -4.0548 3.0216 -1.342 0.180437
## M4M2_M4 -2.6232 3.0401 -0.863 0.388783
## M4M2_M5 -2.0649 3.1557 -0.654 0.513314
## M4M2_M6 -1.1986 3.2377 -0.370 0.711438
## M4M3_M3 -1.9095 3.4804 -0.549 0.583579
## M4M3_M4 -3.3161 3.1422 -1.055 0.291958
## M4M3_M5 -0.9641 3.6624 -0.263 0.792519
## M4M3_M6 -5.1308 3.7245 -1.378 0.169169
## M4M4_M4 0.7499 3.4296 0.219 0.827031
## M4M4_M5 -1.7700 3.6168 -0.489 0.624854
## M4M4_M6 -6.6720 3.7609 -1.774 0.076879 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.214 on 369 degrees of freedom
## Multiple R-squared: 0.9021, Adjusted R-squared: 0.8949
## F-statistic: 125.9 on 27 and 369 DF, p-value: < 2.2e-16
summary(lm(weight~sex+sire+M5,data=markers))
##
## Call:
## lm(formula = weight ~ sex + sire + M5, data = markers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1875 -3.1500 -0.0774 3.5003 13.5327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 305.7865 1.8057 169.347 < 2e-16 ***
## sexM 46.1594 0.5183 89.065 < 2e-16 ***
## siresire10 14.8838 1.1756 12.661 < 2e-16 ***
## siresire2 10.6582 1.2921 8.249 2.86e-15 ***
## siresire3 18.1316 1.2657 14.325 < 2e-16 ***
## siresire4 6.2499 1.1896 5.254 2.53e-07 ***
## siresire5 13.3783 1.1793 11.344 < 2e-16 ***
## siresire6 9.6261 1.2760 7.544 3.58e-13 ***
## siresire7 13.0775 1.1806 11.077 < 2e-16 ***
## siresire8 5.6781 1.3267 4.280 2.39e-05 ***
## siresire9 1.3659 1.1822 1.155 0.2487
## M5M1_M2 1.5100 1.8018 0.838 0.4025
## M5M1_M3 2.3341 1.6972 1.375 0.1699
## M5M1_M4 3.3998 1.7923 1.897 0.0586 .
## M5M1_M5 1.9259 1.9212 1.002 0.3168
## M5M1_M6 2.2328 1.8499 1.207 0.2282
## M5M2_M2 -16.2519 2.1204 -7.665 1.60e-13 ***
## M5M2_M3 -14.5825 1.7261 -8.448 6.90e-16 ***
## M5M2_M4 -14.0451 1.7926 -7.835 5.06e-14 ***
## M5M2_M5 -14.2710 1.9344 -7.378 1.07e-12 ***
## M5M2_M6 -16.9865 2.0890 -8.131 6.51e-15 ***
## M5M3_M3 -14.1392 1.8677 -7.571 3.01e-13 ***
## M5M3_M4 -13.5688 1.6838 -8.058 1.09e-14 ***
## M5M3_M5 -12.4149 1.8419 -6.740 6.10e-11 ***
## M5M3_M6 -14.0551 1.7862 -7.869 4.01e-14 ***
## M5M4_M4 -9.6016 2.0063 -4.786 2.47e-06 ***
## M5M4_M5 -12.3662 1.9930 -6.205 1.47e-09 ***
## M5M4_M6 -10.0045 2.3903 -4.185 3.56e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.04 on 369 degrees of freedom
## Multiple R-squared: 0.9631, Adjusted R-squared: 0.9604
## F-statistic: 357 on 27 and 369 DF, p-value: < 2.2e-16
plot.design(weight~sex+sire+M5,data=markers,col= "red")
anova(lm(weight~sex+sire+M1,data=markers))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 205935 205935 2987.1633 <2e-16 ***
## sire 9 21979 2442 35.4243 <2e-16 ***
## M1 18 958 53 0.7722 0.7331
## Residuals 368 25370 69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(weight~sex+sire+M2,data=markers))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 205935 205935 3008.1491 <2e-16 ***
## sire 9 21979 2442 35.6732 <2e-16 ***
## M2 17 1067 63 0.9166 0.5545
## Residuals 369 25261 68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(weight~sex+sire+M3,data=markers))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 205935 205935 2995.3232 <2e-16 ***
## sire 9 21979 2442 35.5211 <2e-16 ***
## M3 17 959 56 0.8202 0.6698
## Residuals 369 25370 69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(weight~sex+sire+M4,data=markers))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 205935 205935 3051.9725 <2e-16 ***
## sire 9 21979 2442 36.1929 <2e-16 ***
## M4 17 1430 84 1.2462 0.2256
## Residuals 369 24899 67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(weight~sex+sire+M5,data=markers))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 205935 205935 8106.770 < 2.2e-16 ***
## sire 9 21979 2442 96.137 < 2.2e-16 ***
## M5 17 16954 997 39.260 < 2.2e-16 ***
## Residuals 369 9374 25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(car)
Anova(lm(weight~sex+sire+M1+M2+M3+M4+M5, data=markers), type=3)
## Anova Table (Type III tests)
##
## Response: weight
## Sum Sq Df F value Pr(>F)
## (Intercept) 143615 1 5574.5479 <2e-16 ***
## sex 162947 1 6324.9082 <2e-16 ***
## sire 7529 9 32.4717 <2e-16 ***
## M1 351 18 0.7560 0.7510
## M2 293 17 0.6699 0.8320
## M3 347 17 0.7913 0.7031
## M4 574 17 1.3095 0.1845
## M5 14151 17 32.3113 <2e-16 ***
## Residuals 7729 300
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(nlme)
linear=coefficients(lm(weight~sex+sire+M5, data=markers))[c(1,2,12:28)]
random=fixef(lme(weight~sex+M5,random=~1|sire, data=markers))
linear=data.frame(effect=names(linear),fixed=linear)
random=data.frame(effect=names(random),random=random)
comparison=merge(linear,random,by="effect")
comparison=data.frame(comparison, difference=comparison$fixed-comparison$random)
print(comparison)
## effect fixed random difference
## 1 (Intercept) 305.786519 315.102970 -9.316451527
## 2 M5M1_M2 1.510027 1.496079 0.013948130
## 3 M5M1_M3 2.334150 2.320083 0.014067092
## 4 M5M1_M4 3.399809 3.447528 -0.047719235
## 5 M5M1_M5 1.925917 1.968932 -0.043015397
## 6 M5M1_M6 2.232801 2.254522 -0.021720830
## 7 M5M2_M2 -16.251889 -16.350045 0.098156469
## 8 M5M2_M3 -14.582450 -14.609071 0.026620339
## 9 M5M2_M4 -14.045079 -14.090447 0.045368309
## 10 M5M2_M5 -14.271022 -14.315654 0.044631473
## 11 M5M2_M6 -16.986481 -17.019542 0.033060629
## 12 M5M3_M3 -14.139227 -14.165658 0.026430752
## 13 M5M3_M4 -13.568813 -13.597116 0.028302729
## 14 M5M3_M5 -12.414888 -12.437924 0.023035283
## 15 M5M3_M6 -14.055094 -14.085797 0.030703298
## 16 M5M4_M4 -9.601594 -9.611857 0.010263276
## 17 M5M4_M5 -12.366240 -12.302824 -0.063415133
## 18 M5M4_M6 -10.004451 -10.063377 0.058925858
## 19 sexM 46.159432 46.168063 -0.008631174