Chickweight DOE descrpitive Stat
Mudry Stater
2025-02-26
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(lattice)
library(car)## Le chargement a nécessité le package : carDatalibrary(gplots)## Warning: le package 'gplots' a été compilé avec la version R 4.2.3##
## Attachement du package : 'gplots'## L'objet suivant est masqué depuis 'package:stats':
##
## lowesslibrary(psych)##
## Attachement du package : 'psych'## L'objet suivant est masqué depuis 'package:car':
##
## logitlibrary(nlme)
library(nlmeU)## Warning: le package 'nlmeU' a été compilé avec la version R 4.2.3##
## Attachement du package : 'nlmeU'## L'objet suivant est masqué depuis 'package:stats':
##
## sigmaplot(ChickWeight)
str(ChickWeight)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 578 obs. of 4 variables:
## $ weight: num 42 51 59 64 76 93 106 125 149 171 ...
## $ Time : num 0 2 4 6 8 10 12 14 16 18 ...
## $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Diet : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "outer")=Class 'formula' language ~Diet
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "labels")=List of 2
## ..$ x: chr "Time"
## ..$ y: chr "Body weight"
## - attr(*, "units")=List of 2
## ..$ x: chr "(days)"
## ..$ y: chr "(gm)"C=ChickWeight
plot(C$weight~C$Time,col=as.integer(C$Chick))
boxplot(C$weight~C$Time)
###
library(lattice)
##DOE
##according to https://books.google.ch/books?id=XsGX6Jgzo-IC&pg=PA75&hl=fr&source=gbs_toc_r&cad=2#v=onepage&q&f=false DIET control DIET1 10% protein 2%20 % prot 40%prot
#data repèorting deocnt telle us much: as in ex.does ths feediing oof protein is in replacment or added as calories to the food? Unclera from researcher report (Said replacment so we must assume that the total calories or food waeight /daily is the same but i not precising it make some unclear stament and hence conclusion0
#Document your data Know your Data Know your experiemnts)
#TRAJECTORIES GRWOTH CURVES?
xyplot(C$weight~jitter(C$Time)|C$Diet,type=c("p","r","smooth"),pch=20,lwd=2,aspect = 1.53)##Only group 3 depart from linearity in time but mild
xyplot(C$weight~factor(C$Time),groups = C$Chick,type="l")
xyplot(C$weight~factor(C$Time)|C$Diet,groups = C$Chick,type="l")
sum(is.na(C))#Cautious in long fprmat NA can be oomiited by row recoding/removing## [1] 0##chick by diet balanced DOE?
xtabs(~C$Diet)## C$Diet
## 1 2 3 4
## 220 120 120 118countdiet=NULL
#how your code can be wrong
countdiet## NULLunique(C$Chick[C$Diet==1])## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48unique(C$Chick[C$Diet==2])## [1] 21 22 23 24 25 26 27 28 29 30
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48unique(C$Chick[C$Diet==3])## [1] 31 32 33 34 35 36 37 38 39 40
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48unique(C$Chick[C$Diet==4])## [1] 41 42 43 44 45 46 47 48 49 50
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48C$Chick## [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3
## [26] 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5
## [51] 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7
## [76] 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9
## [101] 9 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10 10 10 10 11 11 11 11 11 11
## [126] 11 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13
## [151] 13 13 13 13 13 14 14 14 14 14 14 14 14 14 14 14 14 15 15 15 15 15 15 15 15
## [176] 16 16 16 16 16 16 16 17 17 17 17 17 17 17 17 17 17 17 17 18 18 19 19 19 19
## [201] 19 19 19 19 19 19 19 19 20 20 20 20 20 20 20 20 20 20 20 20 21 21 21 21 21
## [226] 21 21 21 21 21 21 21 22 22 22 22 22 22 22 22 22 22 22 22 23 23 23 23 23 23
## [251] 23 23 23 23 23 23 24 24 24 24 24 24 24 24 24 24 24 24 25 25 25 25 25 25 25
## [276] 25 25 25 25 25 26 26 26 26 26 26 26 26 26 26 26 26 27 27 27 27 27 27 27 27
## [301] 27 27 27 27 28 28 28 28 28 28 28 28 28 28 28 28 29 29 29 29 29 29 29 29 29
## [326] 29 29 29 30 30 30 30 30 30 30 30 30 30 30 30 31 31 31 31 31 31 31 31 31 31
## [351] 31 31 32 32 32 32 32 32 32 32 32 32 32 32 33 33 33 33 33 33 33 33 33 33 33
## [376] 33 34 34 34 34 34 34 34 34 34 34 34 34 35 35 35 35 35 35 35 35 35 35 35 35
## [401] 36 36 36 36 36 36 36 36 36 36 36 36 37 37 37 37 37 37 37 37 37 37 37 37 38
## [426] 38 38 38 38 38 38 38 38 38 38 38 39 39 39 39 39 39 39 39 39 39 39 39 40 40
## [451] 40 40 40 40 40 40 40 40 40 40 41 41 41 41 41 41 41 41 41 41 41 41 42 42 42
## [476] 42 42 42 42 42 42 42 42 42 43 43 43 43 43 43 43 43 43 43 43 43 44 44 44 44
## [501] 44 44 44 44 44 44 45 45 45 45 45 45 45 45 45 45 45 45 46 46 46 46 46 46 46
## [526] 46 46 46 46 46 47 47 47 47 47 47 47 47 47 47 47 47 48 48 48 48 48 48 48 48
## [551] 48 48 48 48 49 49 49 49 49 49 49 49 49 49 49 49 50 50 50 50 50 50 50 50 50
## [576] 50 50 50
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48summary(C$Chick)##note that CHICK 18,16,15,8,44 have incomplete repated measure## 18 16 15 13 9 20 10 8 17 19 4 6 11 3 1 12 2 5 14 7 24 30 22 23 27 28
## 2 7 8 12 12 12 12 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## 26 25 29 21 33 37 36 31 39 38 32 40 34 35 44 45 43 41 47 49 46 50 42 48
## 12 12 12 12 12 12 12 12 12 12 12 12 12 12 10 12 12 12 12 12 12 12 12 12#DOIET 1 control = 20 count
#DIET2,3,4=10
#5' chicks
#OUTLIERS IN DIET2 No weight increase
C$Chick[C$weight<100 & C$Time>10 & C$Diet==2]## [1] 24 24 24 24 24 24
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48plot(C$weight[C$Chick=="24"]~C$Time[C$Chick=="24"],ylim=c(0,250),type="l",add=TRUE,col=2,lwd=2)## Warning in plot.window(...): "add" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "add" n'est pas un paramètre graphique## Warning in axis(side = side, at = at, labels = labels, ...): "add" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "add" n'est pas un
## paramètre graphique## Warning in box(...): "add" n'est pas un paramètre graphique## Warning in title(...): "add" n'est pas un paramètre graphiquepoints(C$weight[C$Chick=="25"]~C$Time[C$Chick=="25"],ylim=c(0,250),type="l")
describe(C$weight[C$Chick=="24"])## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 12 66.25 10.36 70.5 67.7 5.19 42 76 34 -1.15 -0.03 2.99#This one must be removed form df for porper growth curve estimation model
CC=subset(C,C$Chick!=24)
sum((CC$Chick==24)*1)## [1] 0sum(CC$Chick==25)##always check and recheck your code !!!## [1] 12xyplot(CC$weight~factor(CC$Time)|CC$Diet,groups = CC$Chick,type="l")##removed
#Note In group 1 there is also a deth outliers I decided to keep it because we know that in DIET 1 hazard rate seems to be higher despite binomial testing The incertainty hence variability will reflfected by the SE of coeff making contrast with other DIET
#as we dont know what is about DIET regimen supposed to be a controul group???
##There is a SUPER CHICK in Diet 2 (blue) nevertheless the reason I will keep it but for fine tuning modelling growth curve wort to consider removing it!
###efect of scale on graph is not showing what is the truth is:
plot(C$weight[C$Chick=="24"]~C$Time[C$Chick=="24"],type="l",add=TRUE,col=2,lwd=2)## Warning in plot.window(...): "add" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "add" n'est pas un paramètre graphique## Warning in axis(side = side, at = at, labels = labels, ...): "add" n'est pas un
## paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "add" n'est pas un
## paramètre graphique## Warning in box(...): "add" n'est pas un paramètre graphique## Warning in title(...): "add" n'est pas un paramètre graphique
C$Time## [1] 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0
## [26] 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2
## [51] 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4
## [76] 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8
## [101] 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10
## [126] 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12
## [151] 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14
## [176] 0 2 4 6 8 10 12 0 2 4 6 8 10 12 14 16 18 20 21 0 2 0 2 4 6
## [201] 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8
## [226] 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10
## [251] 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12
## [276] 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14
## [301] 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16
## [326] 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18
## [351] 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20
## [376] 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21
## [401] 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0
## [426] 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2
## [451] 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4
## [476] 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6
## [501] 8 10 12 14 16 18 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12
## [526] 14 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14
## [551] 16 18 20 21 0 2 4 6 8 10 12 14 16 18 20 21 0 2 4 6 8 10 12 14 16
## [576] 18 20 21C$Timefac=factor(C$Time)
plot(C$weight~C$Time)
monaif=lm(C$weight~C$Time)
abline(lm(C$weight~C$Time),col=2)
qqnorm(monaif$residuals)
qqline(monaif$residuals)
plot(lag(monaif$residuals),type="h")
lag(monaif$residuals)## 1 2 3 4 5
## 14.53257485 5.92649631 -3.67958222 -16.28566076 -21.89173929
## 6 7 8 9 10
## -22.49781783 -27.10389636 -25.70997490 -19.31605343 -14.92213197
## 11 12 13 14 15
## -4.52821050 -7.33124977 12.53257485 3.92649631 -4.67958222
## 16 17 18 19 20
## -8.28566076 -13.89173929 -12.49781783 -11.10389636 -12.70997490
## 21 22 23 24 25
## -6.31605343 1.07786803 5.47178950 2.66875023 15.53257485
## 26 27 28 29 30
## -6.07350369 -7.67958222 -13.28566076 -13.89173929 -16.49781783
## 31 32 33 34 35
## -18.10389636 -12.70997490 -5.31605343 1.07786803 -5.52821050
## 36 37 38 39 40
## -10.33124977 14.53257485 3.92649631 -6.67958222 -13.28566076
## 41 42 43 44 45
## -23.89173929 -28.49781783 -31.10389636 -42.70997490 -32.31605343
## 46 47 48 49 50
## -31.92213197 -43.52821050 -55.33124977 13.53257485 -3.07350369
## 51 52 53 54 55
## -14.67958222 -20.28566076 -18.89173929 -9.49781783 7.89610364
## 56 57 58 59 60
## 13.29002510 28.68394657 13.07786803 16.47178950 10.66875023
## 61 62 63 64 65
## 13.53257485 3.92649631 -3.67958222 -6.28566076 -0.89173929
## 66 67 68 69 70
## 8.50218217 7.89610364 -2.70997490 -13.31605343 -25.92213197
## 71 72 73 74 75
## -43.52821050 -55.33124977 13.53257485 3.92649631 -5.67958222
## 76 77 78 79 80
## -9.28566076 -8.89173929 -3.49781783 12.89610364 23.29002510
## 81 82 83 84 85
## 49.68394657 64.07786803 84.47178950 92.66875023 14.53257485
## 86 87 88 89 90
## 4.92649631 -1.67958222 -9.28566076 -13.89173929 -22.49781783
## 91 92 93 94 95
## -23.10389636 -34.70997490 -42.31605343 -51.92213197 -78.52821050
## 96 97 98 99 100
## 14.53257485 5.92649631 -3.67958222 -12.28566076 -12.89173929
## 101 102 103 104 105
## -19.49781783 -43.10389636 -58.70997490 -75.31605343 -85.92213197
## 106 107 108 109 110
## -103.52821050 -114.33124977 13.53257485 -1.07350369 -10.67958222
## 111 112 113 114 115
## -17.28566076 -23.89173929 -34.49781783 -44.10389636 -54.70997490
## 116 117 118 119 120
## -67.31605343 -73.92213197 -83.52821050 -88.33124977 15.53257485
## 121 122 123 124 125
## 5.92649631 0.32041778 3.71433924 14.10826071 23.50218217
## 126 127 128 129 130
## 34.89610364 26.29002510 13.68394657 -1.92213197 -22.52821050
## 131 132 133 134 135
## -37.33124977 13.53257485 3.92649631 -6.67958222 -18.28566076
## 136 137 138 139 140
## -25.89173929 -27.49781783 -14.10389636 -15.70997490 -6.31605343
## 141 142 143 144 145
## -0.92213197 -8.52821050 -7.33124977 13.53257485 2.92649631
## 146 147 148 149 150
## -9.67958222 -20.28566076 -32.89173929 -48.49781783 -62.10389636
## 151 152 153 154 155
## -80.70997490 -97.31605343 -104.92213197 -112.52821050 -116.33124977
## 156 157 158 159 160
## 13.53257485 3.92649631 -0.67958222 -1.28566076 3.10826071
## 161 162 163 164 165
## 12.50218217 30.89610364 41.29002510 58.68394657 62.07786803
## 166 167 168 169 170
## 55.47178950 53.66875023 13.53257485 3.92649631 -6.67958222
## 171 172 173 174 175
## -16.28566076 -29.89173929 -47.49781783 -66.10389636 -82.70997490
## 176 177 178 179 180
## 13.53257485 -0.07350369 -13.67958222 -29.28566076 -40.89173929
## 181 182 183 184 185
## -64.49781783 -79.10389636 14.53257485 5.92649631 -1.67958222
## 186 187 188 189 190
## -8.28566076 -14.89173929 -26.49781783 -35.10389636 -47.70997490
## 191 192 193 194 195
## -55.31605343 -62.92213197 -70.52821050 -70.33124977 11.53257485
## 196 197 198 199 200
## -10.07350369 15.53257485 2.92649631 -7.67958222 -18.28566076
## 201 202 203 204 205
## -32.89173929 -44.49781783 -51.10389636 -62.70997490 -62.31605343
## 206 207 208 209 210
## -65.92213197 -59.52821050 -55.33124977 13.53257485 1.92649631
## 211 212 213 214 215
## -8.67958222 -22.28566076 -32.89173929 -42.49781783 -56.10389636
## 216 217 218 219 220
## -61.70997490 -70.31605343 -78.92213197 -88.52821050 -95.33124977
## 221 222 223 224 225
## 12.53257485 4.92649631 -0.67958222 5.71433924 27.10826071
## 226 227 228 229 230
## 47.50218217 83.89610364 89.29002510 106.68394657 121.07786803
## 231 232 233 234 235
## 114.47178950 118.66875023 13.53257485 9.92649631 1.32041778
## 236 237 238 239 240
## -3.28566076 -7.89173929 -20.49781783 -25.10389636 -39.70997490
## 241 242 243 244 245
## -37.31605343 -37.92213197 -39.52821050 -45.33124977 15.53257485
## 246 247 248 249 250
## 6.92649631 -1.67958222 -7.28566076 -7.89173929 -12.49781783
## 251 252 253 254 255
## -6.10389636 -15.70997490 -23.31605343 -22.92213197 -33.52821050
## 256 257 258 259 260
## -37.33124977 14.53257485 6.92649631 -4.67958222 -6.28566076
## 261 262 263 264 265
## -31.89173929 -47.49781783 -63.10389636 -79.70997490 -96.31605343
## 266 267 268 269 270
## -113.92213197 -127.52821050 -138.33124977 12.53257485 3.92649631
## 271 272 273 274 275
## -0.67958222 -2.28566076 4.10826071 8.50218217 12.89610364
## 276 277 278 279 280
## 13.29002510 28.68394657 45.07786803 55.47178950 52.66875023
## 281 282 283 284 285
## 14.53257485 2.92649631 -5.67958222 -6.28566076 -4.89173929
## 286 287 288 289 290
## -1.49781783 2.89610364 -3.70997490 0.68394657 19.07786803
## 291 292 293 294 295
## 32.47178950 38.66875023 11.53257485 0.92649631 -4.67958222
## 296 297 298 299 300
## -7.28566076 -10.89173929 -15.49781783 -18.10389636 -27.70997490
## 301 302 303 304 305
## -24.31605343 -22.92213197 -18.52821050 -20.33124977 11.53257485
## 306 307 308 309 310
## 0.92649631 -4.67958222 -7.28566076 -5.89173929 -1.49781783
## 311 312 313 314 315
## 11.89610364 5.29002510 15.68394657 21.07786803 8.47178950
## 316 317 318 319 320
## 20.66875023 11.53257485 2.92649631 -3.67958222 -6.28566076
## 321 322 323 324 325
## -10.89173929 -9.49781783 0.89610364 -0.70997490 18.68394657
## 326 327 328 329 330
## 44.07786803 75.47178950 96.66875023 14.53257485 2.92649631
## 331 332 333 334 335
## -3.67958222 -8.28566076 -12.89173929 -17.49781783 -18.10389636
## 336 337 338 339 340
## -28.70997490 -25.31605343 -34.92213197 -46.52821050 -62.33124977
## 341 342 343 344 345
## 14.53257485 7.92649631 -0.67958222 -7.28566076 -12.89173929
## 346 347 348 349 350
## -13.49781783 -10.10389636 -12.70997490 1.68394657 18.07786803
## 351 352 353 354 355
## 31.47178950 43.66875023 13.53257485 3.92649631 2.32041778
## 356 357 358 359 360
## 1.71433924 9.10826071 13.50218217 25.89610364 28.29002510
## 361 362 363 364 365
## 52.68394657 77.07786803 87.47178950 92.66875023 11.53257485
## 366 367 368 369 370
## 4.92649631 0.32041778 -3.28566076 -1.89173929 -4.49781783
## 371 372 373 374 375
## 3.89610364 -6.70997490 -17.31605343 -39.92213197 -47.52821050
## 376 377 378 379 380
## -65.33124977 13.53257485 3.92649631 0.32041778 4.71433924
## 381 382 383 384 385
## 9.10826071 18.50218217 30.89610364 35.29002510 66.68394657
## 386 387 388 389 390
## 108.07786803 123.47178950 128.66875023 13.53257485 7.92649631
## 391 392 393 394 395
## 1.32041778 6.71433924 25.10826071 42.50218217 67.89610364
## 396 397 398 399 400
## 87.29002510 118.68394657 146.07786803 157.47178950 160.66875023
## 401 402 403 404 405
## 11.53257485 2.92649631 -1.67958222 -4.28566076 0.10826071
## 406 407 408 409 410
## 0.50218217 11.89610364 15.29002510 29.68394657 41.07786803
## 411 412 413 414 415
## 21.47178950 7.66875023 13.53257485 2.92649631 -6.67958222
## 416 417 418 419 420
## -12.28566076 -17.89173929 -32.49781783 -30.10389636 -38.70997490
## 421 422 423 424 425
## -33.31605343 -28.92213197 -34.52821050 -34.33124977 13.53257485
## 426 427 428 429 430
## 3.92649631 -1.67958222 -6.28566076 0.10826071 -6.49781783
## 431 432 433 434 435
## -5.10389636 3.29002510 23.68394657 46.07786803 76.47178950
## 436 437 438 439 440
## 77.66875023 14.53257485 4.92649631 -1.67958222 -2.28566076
## 441 442 443 444 445
## -8.89173929 -6.49781783 -3.10389636 -4.70997490 1.68394657
## 446 447 448 449 450
## 28.07786803 46.47178950 59.66875023 13.53257485 9.92649631
## 451 452 453 454 455
## 3.32041778 -1.28566076 3.10826071 4.50218217 20.89610364
## 456 457 458 459 460
## 31.29002510 46.68394657 76.07786803 91.47178950 108.66875023
## 461 462 463 464 465
## 14.53257485 5.92649631 3.32041778 4.71433924 5.10826071
## 466 467 468 469 470
## 8.50218217 21.89610364 2.29002510 6.68394657 -1.92213197
## 471 472 473 474 475
## -4.52821050 -8.33124977 14.53257485 3.92649631 0.32041778
## 476 477 478 479 480
## 3.71433924 5.10826071 10.50218217 26.89610364 23.29002510
## 481 482 483 484 485
## 35.68394657 48.07786803 65.47178950 68.66875023 14.53257485
## 486 487 488 489 490
## 9.92649631 6.32041778 15.71433924 33.10826071 41.50218217
## 491 492 493 494 495
## 50.89610364 37.29002510 28.68394657 12.07786803 -4.52821050
## 496 497 498 499 500
## -12.33124977 14.53257485 5.92649631 2.32041778 5.71433924
## 501 502 503 504 505
## 5.10826071 2.50218217 -6.10389636 -12.70997490 -23.31605343
## 506 507 508 509 510
## -39.92213197 13.53257485 4.92649631 -1.67958222 -2.28566076
## 511 512 513 514 515
## 0.10826071 1.50218217 1.89610364 -9.70997490 -21.31605343
## 516 517 518 519 520
## -11.92213197 -6.52821050 -16.33124977 12.53257485 6.92649631
## 521 522 523 524 525
## -0.67958222 1.71433924 3.10826071 4.50218217 10.89610364
## 526 527 528 529 530
## 5.29002510 4.68394657 24.07786803 27.47178950 25.66875023
## 531 532 533 534 535
## 13.53257485 7.92649631 3.32041778 -1.28566076 2.10826071
## 536 537 538 539 540
## 7.50218217 14.89610364 6.29002510 -0.31605343 -0.92213197
## 541 542 543 544 545
## 6.47178950 -7.33124977 11.53257485 4.92649631 -0.67958222
## 546 547 548 549 550
## -0.28566076 6.10826071 9.50218217 20.89610364 19.29002510
## 551 552 553 554 555
## 53.68394657 75.07786803 99.47178950 109.66875023 12.53257485
## 556 557 558 559 560
## 7.92649631 1.32041778 4.71433924 10.10826071 12.50218217
## 561 562 563 564 565
## 18.89610364 15.29002510 15.68394657 17.07786803 29.47178950
## 566 567 568 569 570
## 24.66875023 13.53257485 8.92649631 4.32041778 3.71433924
## 571 572 573 574 575
## 7.10826071 6.50218217 21.89610364 24.29002510 36.68394657
## 576 577 578
## 48.07786803 60.47178950 51.66875023
## attr(,"tsp")
## [1] 0 577 1plot(monaif)
barplot(table(C$Chick,C$Timefac),main="Count of mesures by time",xlab="Time in days",ylab="daily measure counts")##count of reepeted mesures 
#mortaility is nt reported cause of death ,ight be natural check in gorups but no firm comclusion witrh low n by treatment
missing=tapply(C$weight,list(C$Time,C$Diet),FUN=function(x) length(x[!is.na(x)]))
missing=as.data.frame(missing)
par(mfrow=c(2,2))
plot(missing$`1`,type="s",col=2,main="Chick death loss DIET1",ylim=c(8,20))
plot(missing$`2`,type="s",col=2,main="Chick death loss DIET2",ylim=c(8,20))
plot(missing$`3`,type="s",col=2,main="Chick death loss DIET3",ylim=c(8,20))
plot(missing$`4`,type="s",col=2,main="Chick death loss DIET4",ylim=c(8,20))
lossprob=5/length(unique(C$Chick))
lossprob## [1] 0.1binom.test(1,12,p=0.1) ###loss on DIET 4 is not significant from random chance mother nature##
## Exact binomial test
##
## data: 1 and 12
## number of successes = 1, number of trials = 12, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.1
## 95 percent confidence interval:
## 0.002107593 0.384796165
## sample estimates:
## probability of success
## 0.08333333#as conclusion DIET1 cause hazard death death rate than random chance throughout stratum of DIET
binom.test(4,20,p=0.1,alternative="greater") ###loss on DIET 4 is not significant from random chance mother nature#here we are at p= 0.20 stratum 1 loss;##
## Exact binomial test
##
## data: 4 and 20
## number of successes = 4, number of trials = 20, p-value = 0.133
## alternative hypothesis: true probability of success is greater than 0.1
## 95 percent confidence interval:
## 0.07135388 1.00000000
## sample estimates:
## probability of success
## 0.2#Although the size of population has doubled the p hat has doubled too wich is a signed of a possible changes un hazard rate although the p value is nit signif but has lowered
#RECOMMENDATION :INVESTIGATION of cause of Death must be done on DIET1
#sphericity
describe(C$weight)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 578 121.82 71.07 103 113.18 69.68 35 373 338 0.96 0.34 2.96#range from day 0 to 21
338/35## [1] 9.657143# 9,7 his weight
dotchart(tapply(C$weight,C$Chick,sd),groups=as.numeric(C$Diet),cex=0.8,gcolor=rainbow(4))## Warning in dotchart(tapply(C$weight, C$Chick, sd), groups =
## as.numeric(C$Diet), : 'x' n'est ni un vecteur ni une matrice : utilisation de
## as.numeric(x)## Warning in seq_len(n) + 2 * of.1: la taille d'un objet plus long n'est pas
## multiple de la taille d'un objet plus courtabline(v=71,lty=3,lwd=0.7)
C1=split(C,C$Diet,drop=F)
str(C1)## List of 4
## $ 1:Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 220 obs. of 5 variables:
## ..$ weight : num [1:220] 42 51 59 64 76 93 106 125 149 171 ...
## ..$ Time : num [1:220] 0 2 4 6 8 10 12 14 16 18 ...
## ..$ Chick : Ord.factor w/ 20 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
## ..$ Diet : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ Timefac: Factor w/ 12 levels "0","2","4","6",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..- attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "labels")=List of 2
## .. ..$ x: chr "Time"
## .. ..$ y: chr "Body weight"
## ..- attr(*, "units")=List of 2
## .. ..$ x: chr "(days)"
## .. ..$ y: chr "(gm)"
## ..- attr(*, "outer")=Class 'formula' language ~Diet
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "FUN")=function (x)
## ..- attr(*, "order.groups")= logi TRUE
## $ 2:Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 120 obs. of 5 variables:
## ..$ weight : num [1:120] 40 50 62 86 125 163 217 240 275 307 ...
## ..$ Time : num [1:120] 0 2 4 6 8 10 12 14 16 18 ...
## ..$ Chick : Ord.factor w/ 10 levels "24"<"30"<"22"<..: 10 10 10 10 10 10 10 10 10 10 ...
## ..$ Diet : Factor w/ 1 level "2": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ Timefac: Factor w/ 12 levels "0","2","4","6",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..- attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "labels")=List of 2
## .. ..$ x: chr "Time"
## .. ..$ y: chr "Body weight"
## ..- attr(*, "units")=List of 2
## .. ..$ x: chr "(days)"
## .. ..$ y: chr "(gm)"
## ..- attr(*, "outer")=Class 'formula' language ~Diet
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "FUN")=function (x)
## ..- attr(*, "order.groups")= logi TRUE
## $ 3:Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 120 obs. of 5 variables:
## ..$ weight : num [1:120] 42 53 62 73 85 102 123 138 170 204 ...
## ..$ Time : num [1:120] 0 2 4 6 8 10 12 14 16 18 ...
## ..$ Chick : Ord.factor w/ 10 levels "33"<"37"<"36"<..: 4 4 4 4 4 4 4 4 4 4 ...
## ..$ Diet : Factor w/ 1 level "3": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ Timefac: Factor w/ 12 levels "0","2","4","6",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..- attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "labels")=List of 2
## .. ..$ x: chr "Time"
## .. ..$ y: chr "Body weight"
## ..- attr(*, "units")=List of 2
## .. ..$ x: chr "(days)"
## .. ..$ y: chr "(gm)"
## ..- attr(*, "outer")=Class 'formula' language ~Diet
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "FUN")=function (x)
## ..- attr(*, "order.groups")= logi TRUE
## $ 4:Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 118 obs. of 5 variables:
## ..$ weight : num [1:118] 42 51 66 85 103 124 155 153 175 184 ...
## ..$ Time : num [1:118] 0 2 4 6 8 10 12 14 16 18 ...
## ..$ Chick : Ord.factor w/ 10 levels "44"<"45"<"43"<..: 4 4 4 4 4 4 4 4 4 4 ...
## ..$ Diet : Factor w/ 1 level "4": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ Timefac: Factor w/ 12 levels "0","2","4","6",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..- attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "labels")=List of 2
## .. ..$ x: chr "Time"
## .. ..$ y: chr "Body weight"
## ..- attr(*, "units")=List of 2
## .. ..$ x: chr "(days)"
## .. ..$ y: chr "(gm)"
## ..- attr(*, "outer")=Class 'formula' language ~Diet
## .. .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## ..- attr(*, "FUN")=function (x)
## ..- attr(*, "order.groups")= logi TRUEGroup1=C1$`1`
Group2=C1$`2`
Group3=C1$`3`
Group4=C1$`4`
par(mfrow=c(2,2))
plot(tapply(Group1$weight,Group1$Chick,sd),ylim=c(0,140),type="h")
plot(tapply(Group2$weight,Group2$Chick,sd),ylim=c(0,140),type="h")
plot(tapply(Group3$weight,Group3$Chick,sd),ylim=c(0,140),type="h")
plot(tapply(Group4$weight,Group4$Chick,sd),ylim=c(0,140),type="h")
dev.off()## null device
## 1plot(tapply(C$weight,C$Diet,sd),type="h")
plot(tapply(C$weight,C$Diet,sd),type="h",ylim=c(0,140))
abline(h=71.2,lty=3,lwd=0.7,col=3)
#Empirical Within-Subject Residuals
C$weight[C$Chick==18]## [1] 39 35C$weight## [1] 42 51 59 64 76 93 106 125 149 171 199 205 40 49 58 72 84 103
## [19] 122 138 162 187 209 215 43 39 55 67 84 99 115 138 163 187 198 202
## [37] 42 49 56 67 74 87 102 108 136 154 160 157 41 42 48 60 79 106
## [55] 141 164 197 199 220 223 41 49 59 74 97 124 141 148 155 160 160 157
## [73] 41 49 57 71 89 112 146 174 218 250 288 305 42 50 61 71 84 93
## [91] 110 116 126 134 125 42 51 59 68 85 96 90 92 93 100 100 98 41
## [109] 44 52 63 74 81 89 96 101 112 120 124 43 51 63 84 112 139 168
## [127] 177 182 184 181 175 41 49 56 62 72 88 119 135 162 185 195 205 41
## [145] 48 53 60 65 67 71 70 71 81 91 96 41 49 62 79 101 128 164
## [163] 192 227 248 259 266 41 49 56 64 68 68 67 68 41 45 49 51 57
## [181] 51 54 42 51 61 72 83 89 98 103 113 123 133 142 39 35 43 48
## [199] 55 62 65 71 82 88 106 120 144 157 41 47 54 58 65 73 77 89
## [217] 98 107 115 117 40 50 62 86 125 163 217 240 275 307 318 331 41 55
## [235] 64 77 90 95 108 111 131 148 164 167 43 52 61 73 90 103 127 135
## [253] 145 163 170 175 42 52 58 74 66 68 70 71 72 72 76 74 40 49
## [271] 62 78 102 124 146 164 197 231 259 265 42 48 57 74 93 114 136 147
## [289] 169 205 236 251 39 46 58 73 87 100 115 123 144 163 185 192 39 46
## [307] 58 73 92 114 145 156 184 207 212 233 39 48 59 74 87 106 134 150
## [325] 187 230 279 309 42 48 59 72 85 98 115 122 143 151 157 150 42 53
## [343] 62 73 85 102 123 138 170 204 235 256 41 49 65 82 107 129 159 179
## [361] 221 263 291 305 39 50 63 77 96 111 137 144 151 146 156 147 41 49
## [379] 63 85 107 134 164 186 235 294 327 341 41 53 64 87 123 158 201 238
## [397] 287 332 361 373 39 48 61 76 98 116 145 166 198 227 225 220 41 48
## [415] 56 68 80 83 103 112 135 157 169 178 41 49 61 74 98 109 128 154
## [433] 192 232 280 290 42 50 61 78 89 109 130 146 170 214 250 272 41 55
## [451] 66 79 101 120 154 182 215 262 295 321 42 51 66 85 103 124 155 153
## [469] 175 184 199 204 42 49 63 84 103 126 160 174 204 234 269 281 42 55
## [487] 69 96 131 157 184 188 197 198 199 200 42 51 65 86 103 118 127 138
## [505] 145 146 41 50 61 78 98 117 135 141 147 174 197 196 40 52 62 82
## [523] 101 120 144 156 173 210 231 238 41 53 66 79 100 123 148 157 168 185
## [541] 210 205 39 50 62 80 104 125 154 170 222 261 303 322 40 53 64 85
## [559] 108 128 152 166 184 203 233 237 41 54 67 84 105 122 155 175 205 234
## [577] 264 264xyplot(C$diffwithin~C$Chick)## Warning in diff(as.numeric(y[ord])): NAs introduits lors de la conversion
## automatique## Warning in (function (x, y, type = "p", groups = NULL, pch = if
## (is.null(groups)) plot.symbol$pch else superpose.symbol$pch, : NAs introduits
## lors de la conversion automatiquehead(C)## Grouped Data: weight ~ Time | Chick
## weight Time Chick Diet Timefac
## 1 42 0 1 1 0
## 2 51 2 1 1 2
## 3 59 4 1 1 4
## 4 64 6 1 1 6
## 5 76 8 1 1 8
## 6 93 10 1 1 10xyplot(C$diffwithin~factor(C$Time)|C$Chick,type="l",abline=c(h=0,lwd=2),
col=as.numeric(C$Diet))## Warning in diff(as.numeric(y[ord])): NAs introduits lors de la conversion
## automatique## Warning in (function (x, y, type = "p", groups = NULL, pch = if
## (is.null(groups)) plot.symbol$pch else superpose.symbol$pch, : NAs introduits
## lors de la conversion automatique###probelm some chick has no mean zeto why
summary(C)## weight Time Chick Diet Timefac
## Min. : 35.0 Min. : 0.00 13 : 12 1:220 0 : 50
## 1st Qu.: 63.0 1st Qu.: 4.00 9 : 12 2:120 2 : 50
## Median :103.0 Median :10.00 20 : 12 3:120 4 : 49
## Mean :121.8 Mean :10.72 10 : 12 4:118 6 : 49
## 3rd Qu.:163.8 3rd Qu.:16.00 17 : 12 8 : 49
## Max. :373.0 Max. :21.00 19 : 12 10 : 49
## (Other):506 (Other):282str(C$Chick)## Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...unique(C$Chick )## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48mean(C$weight[C$Chick==18])## [1] 37X=data.frame(C$Diet,C$Chick)
Y=C$weight
plot.design(X,Y,cex=0.8)
interaction.plot(cut(C$Time,2),C$Diet,C$weight,lwd=2,col=c(1,2,3,4),ylim=c(40,230))
interaction.plot(cut(C$Time,8),C$Diet,ylim=c(0,250),C$weight,lwd=2,col=c(1,2,3,4))
interaction.plot(cut(C$Time,8),C$Diet,ylim=c(0,250),C$weight,lwd=2,col=c(1,2,3,4),fun="mean")
##sd by two strata is lower than one factor this is a property as a stratification sampling has lower variance always demo see.Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.