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 only consider the 30 animal data in which only treatment A is given. Karim et al. (2022) already concluded that the intrinsic growth rate r varies with time as r(t) = b exp (−ct)t, where b, c are positive real numbers. Using ISRP computed by maximizing the local likelihood, here we also identify the continuously varying form for the parameter r. We first re-scaling the measurement schedules to t = 1, 2, . . . , 10 and fit the exponential model into the cattle growth data and estimate the ISRP of the growth parameter r. We will get the estimated value of r for 10 − 1 = 9 time points. We plot the ISRP profile of r. It is evident from the ISRP profile that r varies with time. To verify whether it follows r(t) = b exp (−ct)t, we fit the functional form of r into the ISRP profile. This plot clearly indicates that r follows the same varying form of time as Karim et al. (2022).
# Set the working directory to load the real data
setwd("D:/My Paper/2nd Paper MLE vs ISRP/RDA")
cattle_data_1=read.delim("cattle_data.txt",sep=" ")
colnames(cattle_data_1)[1:13]<-c("id","G","0","14","28","42","56","70","84","98","112","126","133")
cattle_data=cattle_data_1[,1:13]
data_B=cattle_data[1:30,2:12] #for Group B
data_A=cattle_data[31:60,2:13] #for Group A
t=c(0,14,28,42,56,70,84,98,112,126,133)
write.table(data_A[,-1], "data_A.txt", row.names = FALSE, col.names = FALSE)
t = seq(from = 1, to = 11, by = 1) # time points
q=length(t)
data=as.matrix(read.table("data_A.txt", header = FALSE))
n= 30
View(data)
# The following parameters will be considered as the fixed set up for simulation.
r= 0.05
x0 = 220
sigma2= 1
rho= 0.1
mu = x0*exp(r*t)
plot(t,mu)
fun_likelihood = function(param){
r = param[1]
sigma2 = param[2]
rho = param[3]
mu = x0*exp(r*t)
cov_mat = matrix(data = NA, nrow = length(t), ncol = length(t))
for (i in 1:q) {
for (j in 1:q) {
cov_mat[i,j] = sigma2*rho^(abs(i-j))
}
}
exponent = 0
for(i in 1:n){
exponent = exponent + t(data[i,]-mu)%*%(solve(cov_mat))%*%(data[i,]-mu)
}
log_lik = - (n*q/2)*log(2*pi) - (n/2)*log(det(cov_mat)) - (1/2)*exponent
return(-log_lik)
}
param_0 = c(0.1, 1, 0.1)
out = optim(param_0, fun_likelihood, method = "L-BFGS-B",
hessian = TRUE, lower = c(0.001, 0.1, -0.001), upper = c(1, 10, 0.999))
out$par
solve(out$hessian)
# Local MLE
est_local = matrix(data = NA, nrow = q-1, ncol = 3)
fun_local_likelihood = function(param){
r = param[1]
sigma2 = param[2]
rho = param[3]
mu = x0*exp(r*t)
cov_mat = matrix(data = NA, nrow = length(t), ncol = length(t))
for (i in 1:q) {
for (j in 1:q) {
cov_mat[i,j] = sigma2*rho^(abs(i-j))
}
}
# localized estimates
mu_local = mu[ind:(ind + 1)]
data_local = data[,ind:(ind + 1)]
cov_mat_local = cov_mat[ind:(ind+1), ind:(ind+1)]
exponent = 0
for(i in 1:n){
exponent = exponent + t(data_local[i,] - mu_local)%*%(solve(cov_mat_local))%*%(data_local[i,] - mu_local)
}
log_lik = - (n*length(ind:(ind+1))/2)*log(2*pi) - (n/2)*log(det(cov_mat_local)) - (1/2)*exponent
return(-log_lik)
}
# likelihood function
for(j in 1:(q-1)){
ind = j
param_0 = c(0.1, 1, 0.1)
out = optim(param_0, fun_local_likelihood, method = "L-BFGS-B",
hessian = TRUE, lower = c(0.001, 0.1, -0.001), upper = c(1, 10, 0.999))
est_local[j, ] = out$par
}
colnames(est_local) = c("r", "sigma2", "rho")
par(mar = c(5.1, 5.1, 4.1, 2.1))
plot(1:(q-1), est_local[,1], type = "b", col = "red", lwd = 2,
ylab = expression(widehat(r(Delta~t))[MLE]),
xlab = "Time")
r_hat=est_local[,1]
x = 1:10
fit = nls(r_hat ~ b*exp(-c*x)*x, start = list(b=1,c=3))
points(x,r_hat, col = "red", lwd = 2, cex= 1.5, pch = 19)
curve(coef(fit)[1]*exp(-coef(fit)[2]*x)*x, lwd = 3, col = "blue", add = TRUE)
legend("bottomright", legend = bquote(r(t) == b*e^(-c*t)*t),
bty = "n", lwd = 3, col = "blue")
b_hat = coef(fit)[1]
c_hat = coef(fit)[2]
AIC(fit)