Cattle Growth Data
Data
REAL Data Set for this analysis:
id G 0 14 28 42 56 70 84 98 112 126 133
1 B 210 215 230 244 259 266 277 292 292 290 264
2 B 230 240 258 277 277 293 300 323 327 340 343
3 B 226 233 248 277 297 313 322 340 354 365 362
4 B 233 239 253 277 292 310 318 333 336 353 338
5 B 238 241 262 282 300 314 319 331 338 348 338
6 B 225 228 237 261 271 288 300 316 319 333 330
7 B 224 225 239 257 268 290 304 313 310 318 318
8 B 237 241 255 276 293 307 312 336 336 344 328
9 B 237 224 234 239 256 266 276 300 302 293 269
10 B 233 239 259 283 294 313 320 347 348 362 352
11 B 217 222 235 256 267 285 295 317 315 308 301
12 B 228 223 246 266 277 287 300 312 308 328 333
13 B 241 247 268 290 309 323 336 348 359 372 370
14 B 221 221 240 253 273 282 292 307 306 317 318
15 B 217 220 235 259 262 276 284 305 303 315 317
16 B 214 221 237 256 271 283 287 314 316 320 298
17 B 224 231 241 256 265 283 295 314 313 328 334
18 B 200 203 221 236 248 262 276 294 291 311 310
19 B 238 232 252 268 285 298 303 320 324 320 327
20 B 230 222 243 253 268 284 290 316 314 330 330
21 B 217 224 242 265 284 302 309 324 328 338 334
22 B 209 209 221 238 256 267 281 295 301 309 289
23 B 224 227 245 267 279 294 312 328 329 297 297
24 B 230 231 244 261 272 283 294 318 320 333 338
25 B 216 218 223 243 259 270 270 290 301 314 297
26 B 231 239 254 276 294 304 317 335 333 319 307
27 B 207 216 228 255 275 285 296 314 319 330 330
28 B 227 236 251 264 276 287 297 315 309 313 294
29 B 221 232 251 274 284 295 300 323 319 333 322
30 B 233 238 254 266 282 294 295 310 320 327 326
31 A 233 224 245 258 271 287 287 287 290 293 297
32 A 231 238 260 273 290 300 311 313 317 321 326
33 A 232 237 245 265 285 298 304 319 317 334 329
34 A 239 246 268 288 308 309 327 324 327 336 341
35 A 215 216 239 264 282 299 307 321 328 332 337
36 A 236 226 242 255 263 277 290 299 300 308 310
37 A 219 229 246 265 279 292 299 299 298 300 290
38 A 231 245 270 292 302 321 322 334 323 337 337
39 A 230 228 243 255 272 276 277 289 289 300 303
40 A 232 240 247 263 275 286 294 302 308 319 326
41 A 234 237 259 289 311 324 342 347 355 368 368
42 A 237 235 258 263 282 304 318 327 336 349 353
43 A 229 234 254 276 294 315 323 341 346 352 357
44 A 220 227 248 273 290 308 322 326 330 342 343
45 A 232 241 255 276 293 309 310 330 326 329 330
46 A 210 225 242 260 272 277 273 295 292 305 306
47 A 229 241 252 265 274 285 303 308 315 328 328
48 A 204 198 217 233 251 258 272 283 279 295 298
49 A 220 221 236 260 274 295 300 301 310 318 316
50 A 233 234 250 268 280 298 308 319 318 336 333
51 A 234 234 254 274 294 306 318 334 343 349 350
52 A 200 207 217 238 252 267 284 282 282 284 288
53 A 220 213 229 252 254 273 293 289 294 292 298
54 A 225 239 254 269 289 308 313 324 327 347 344
55 A 236 245 257 271 294 307 317 327 328 328 325
56 A 231 231 237 261 274 285 291 301 307 315 320
57 A 208 211 238 254 267 287 306 312 320 337 338
58 A 232 248 261 285 292 307 312 323 318 328 329
59 A 233 241 252 273 301 316 332 336 339 348 345
60 A 221 219 231 251 270 272 287 294 292 292 299
Description of the Data
We consider a data set on cattle growth studied by Kenward, 1987. The cattle data contains the weight of individual cattle at eleven-time points over a 133-day period. The animals were given two treatments, namely A, and B for intestinal parasites. We consider the data of 30 animals after giving treatment A. Here we make an attempt to choose the appropriate mean growth profile based on the data. The measurement schedules are rescaled to t = 0 , 1 , 2 , . . . , 10.
Results:
We first plot the ISRP profile of r with time, and the size profile of an individual animal, respectively. Given the average size profile, a sigmoid-shaped growth curve would be an appropriate choice to start the analysis. However, the ISRP profile indicates that the logistic growth model is not appropriate for describing the growth pattern, as the ISRP first increases, and then decreases. So, according to our proposal in this manuscript, we start with the base model as exponential, and the empirical investigation suggests that the growth coefficient r varies as a function of time as r(t) = a exp(−bt) t^(c), where a, b, c are positive. Such a choice of the ISRP of r has been investigated by Bhowmick et al. (2014). Since for the non-integer value of c, the exact solution can not be obtained, we consider the final model for c = 1, and c = 2 only (extended Gompertz model). The logistic model, and the extended Gompertz model ( c = 1, and c = 2 ) have been fitted to the size profiles. The corresponding Akaike Information Criteria (AIC) values for these three models are 68.9423, 46.14059, and 55.42874, respectively. Therefore, the best model turns out to be the exponential model with continuously varying growth coefficients r(t) = a exp(−bt) t. To be more precise with our conclusion, we simulate B = 1000 bootstrap data sets by randomly selecting rows of the data with replacement. All three designated models have been fitted to each bootstrap sample, and corresponding AIC values are recorded. The bootstrap distribution of the AIC for each model gives a strong indication for the selection of the extended Gompertz model with c = 1.
R Code
# select the working directory where you stored the data
setwd("D:/STUDY/PhD_Work/Aktar/Chapter 1")
cattle_data_1=read.delim("cattle_data.txt",sep=" ") # load the data in r for analysis
colnames(cattle_data_1)[1:12]<-c("G","0","14","28","42","56","70","84","98","112","126","133")
cattle_data=cattle_data_1[,1:12]
data_B=cattle_data[1:30,2:12]#for Group B
data_A=cattle_data[31:60,2:12]#for Group B
t=c(0,14,28,42,56,70,84,98,112,126,133)
# Visualization of the data
View(data_A)
View(data_B)
matplot(t(data_A), type = "l")
lines(colMeans(data_A), type = "l", lwd = 3, col = "red")
mean_A = colMeans(data_A)
plot(t, mean_A, col = "red", lwd = 2, pch = 19)
RGR_A = log(mean_A[2:length(mean_A)]/mean_A[1:(length(mean_A)-1)])
plot(t[1:(length(t)-1)], RGR_A, col = "red", lwd = 2, pch = 19)
RGR_A = matrix(data = NA, nrow = nrow(data_A), ncol = (ncol(data_A)-1))
for (i in 1:nrow(data_A)) {
for(j in 1:(ncol(data_A)-1)){
RGR_A[i,j] = log(data_A[i, j+1]/data_A[i,j])
}
}
lines(t[1:(length(t)-1)],colMeans(RGR_A), type = "l", lwd = 3, col = "blue")
# Non linear fitting
y = colMeans(RGR_A); x = 0:(length(y)-1)
matplot(x,t(RGR_A), lwd = 2, col = "grey", type = "l", xlab = "Time (weeks)",
ylab = expression(r[j](Delta~t)), cex.lab = 1.4)
fit = nls(y ~ a*exp(-b*x)*(x^c), start = list(a = 1, b=1,c=3))
points(x,y, col = "red", lwd = 2, cex= 1.5, pch = 19)
curve(coef(fit)[1]*exp(-coef(fit)[2]*x)*x^(coef(fit)[3]),
lwd = 3,
col = "blue", add = TRUE)
legend("topright", legend = bquote(r(t) == a*e^(-b*t)*t^c),
bty = "n", lwd = 3, col = "blue")
a_hat = coef(fit)[1]
b_hat = coef(fit)[2]
c_hat = coef(fit)[3]
AIC(fit)
matplot(0:length(x), t(data_A), type = "l", lwd =2, col="lightgrey",
ylab = "Size", xlab = "Time (weeks)", cex.lab = 1.5)
points(0:length(x),colMeans(data_A), pch = 19, cex = 1.5, lwd =2,
col = "red")
# c= 2 model
t = 0:(length(t)-1)
mean = colMeans(data_A)
fit_log = nls(mean ~ x0/(1+exp(b*(c-t))), start = list(x0=220, b=.34, c=1))
fit_new1 = nls(mean ~ x0*exp((b/c)*(1/c-(exp(-c*t))*(t+1/c))), start = list(x0=220, b=0.1, c=.4))
fit_new2 = nls(mean ~ x0*exp((b/c)*(1/c-(exp(-c*t))*(t^2+2*t/c+2/c^2))), start = list(x0=220, b=0.08, c=0.5))
AIC(fit_log)
AIC(fit_new1)
AIC(fit_new2)
lines(t, fitted.values(fit_new1), col = "blue", lwd = 3, lty = 1 )
lines(t, fitted.values(fit_log), col = "magenta", lwd = 3, lty =2)
lines(t, fitted.values(fit_new2), col = "darkgreen", lwd = 3, lty=3)
legend("topleft", legend = c("Logistic model", "Extended Gompertz (c=2)", "Extended Gompertz (c=1)"),
col = c("magenta", "blue", "darkgreen"), lwd = c(3,3,3), lty = c(2,1, 3), bty = "n")
# Bootstrap analysis
B = 1000
AIC_fit_log = numeric(length = B)
AIC_fit_new1 = numeric(length = B)
AIC_fit_new2 = numeric(length = B)
for (i in 1:B) {
ind = sample(1:nrow(data_A), size = nrow(data_A), replace = TRUE)
data_b = data_A[ind, ]
mean_b = colMeans(data_b)
t = 0:(length(t)-1)
fit_log_b = nls(mean_b ~ x0/(1+exp(b*(c-t))), start = list(x0=220, b=.34, c=1))
fit_new1_b = nls(mean_b ~ x0*exp((b/c)*(1/c-(exp(-c*t))*(t+1/c))), start = list(x0=220, b=0.1, c=.4))
fit_new2_b = nls(mean_b ~ x0*exp((b/c)*(1/c-(exp(-c*t))*(t^2+2*t/c+2/c^2))), start = list(x0=220, b=0.08, c=0.5))
AIC_fit_log[i] = AIC(fit_log_b)
AIC_fit_new1[i] = AIC(fit_new1_b)
AIC_fit_new2[i] = AIC(fit_new2_b)
}
h1 = hist(AIC_fit_log)
h2 = hist(AIC_fit_new1)
h3 = hist(AIC_fit_new2)
par(mfrow=c(1,1))
plot(h1$breaks, c(h1$counts,0), lwd =3, col = "blue", type = "s", xlim = c(35,75), ylim =c(0,350),lty =4, ylab = "Frequency", xlab = "Bootstrap AIC values")
lines(h2$breaks, c(h2$counts, 0), lwd =3 , col = "red", type = "s", lty = 5)
lines(h3$breaks, c(h3$counts, 0), lwd =3 , col = "magenta",
type = "s", lty = 6)
legend("topleft", col = c("blue", "magenta", "red"),
legend = c("Logistic model", "Extended Gompertz (c=2)","Extended Gompertz (c=1)" ),lwd = c(3,3,3), bty = "n", lty = 4:6)