Measuring and Modeling Foster Kitten Growth: R Stats with Cats

cats

Data Updated 2019-07-12

My partner brought home a box of three week old kittens that had been brought to her work. The man who found them said that during a heavy rainstorm heard tiny meows in his yard and then saw a little kitten face scrambling to get out of a depression. He went to rescue the baby kitten and found FIVE kittens, most of them were totally submerged in rising water.

When we got the kittens home, I did the best thing I knew how to do: start collecting data on their growth.

We color coded their tail tips with the nail polish we had on hand which led to their color names.

I made a quick Excel workbook to collect the date, weight in grams, and comments for each cat. A few days later, I forgot to weigh before dinner so I weighed after eating. I realized that it would be a good way to track their food consumption. I’m glad I took the two measures. The kittens got a little older and very wiggly before meals. The measurement error would have another measurement to track average growth.

The First Week of Fostering

The Second Week of Fostering

Now Plot Their Progress

Load some packages and data. I used the forcats package to re-code factors and was looking for a way to use purrr but didn’t have a good use case.

library(knitr)
library(tidyverse)
library(gganimate)
library(sjPlot)

load("Kittens.Rda")

kable(head(kittens))
Date Cat Weight Pre_post.eating Comments
2019-06-08 Blue 215 Pre-eating Goopy eyes. Applied antibiotics.
2019-06-08 Green 249 Pre-eating Gave eye antibiotics
2019-06-08 Lavender 233 Pre-eating NA
2019-06-08 Maroon 252 Pre-eating NA
2019-06-08 No Tip 255 Pre-eating NA
2019-06-09 Green 265 Pre-eating Diarrhea; eye antibiotics

Plotting Kitten Growth

The raw data need a bit of re-coding. I do several things in the next block:

  1. Make a color variable to use in the plots.
  2. Re-code the “weight.taken” variable so that it will be correct in the plots.
    • Order the weight.taken variable levels so that the before eating comes first.
  3. Make a sex variable based on the cat’s sex.
  4. Make a numeric date variable that will be used in the regression. Zero is the first day.
kittens.coded<-kittens %>%
  mutate(color=fct_recode(Cat, 
                          "black" = "No Tip",
                          "violet" = "Lavender" ),
         color=tolower(color),
         weight.taken=fct_recode(Pre_post.eating,
                                 "Before Eating" = "Pre-eating" ,
                                 "After Eating" = "Post-eating"),
         weight.taken=fct_relevel(weight.taken,"Before Eating"),
         sex=fct_collapse(Cat,
                          female= c("Blue", "Maroon"),
                          male= c("Green", "Lavender", "No Tip")),
         day.num=as.numeric(Date)-18055)

Now to plot with ggplot2.

Note: No Tip has some missing days because someone peed on my data collection form before I had a chance to record those days

kittentraj<-kittens.coded %>%
  ggplot(aes(Date, Weight, color=Cat)) +
  geom_line() +
  scale_color_manual(values = kittens.coded$color)+
  facet_grid(~ weight.taken) + 
  ggtitle("Foster Kitten Weight (Grams)") 

kittentraj

It is also helpful to split the plots by each kitten so that we can see their weight pre- and post-eating easier.

kittens.sep<-kittens.coded %>%
  ggplot(aes(Date, Weight, color=weight.taken)) +
  geom_line() +
  facet_wrap(~ Cat) + 
  ggtitle("Foster Kitten Weight (Grams)") + 
  labs(color = "Weight Taken") 

kittens.sep

Animating Plots

I wanted to play around a bit more with the gganimate package so I made gifs of the growth charts.

kittentraj + 
  geom_point(aes(group = seq_along(Date))) +
  transition_reveal(Date)

kittens.sep + 
  geom_point(aes(group = seq_along(Date))) +
  transition_reveal(Date)

Curvilinear Growth Model

In the OLS, I add a polynomial term for day number. The I() function lets you evaluate the term inside the function without the +, *, or ^ being interpreted as a part of the model equation.

weight.curv<-lm(Weight ~ day.num + I(day.num^2) + Pre_post.eating +  Cat, kittens.coded)

summary(weight.curv)
## 
## Call:
## lm(formula = Weight ~ day.num + I(day.num^2) + Pre_post.eating + 
##     Cat, data = kittens.coded)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.977 -17.430   0.809  15.679  86.522 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               238.87616    6.60877  36.145  < 2e-16 ***
## day.num                    12.50165    0.70128  17.827  < 2e-16 ***
## I(day.num^2)                0.21268    0.01971  10.788  < 2e-16 ***
## Pre_post.eatingPre-eating -22.71086    3.34772  -6.784 8.07e-11 ***
## CatGreen                   57.25390    5.17651  11.060  < 2e-16 ***
## CatLavender                53.15350    5.17651  10.268  < 2e-16 ***
## CatMaroon                  20.38889    5.17614   3.939 0.000106 ***
## CatNo Tip                  45.13954    5.34632   8.443 2.33e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.9 on 256 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.9809, Adjusted R-squared:  0.9804 
## F-statistic:  1878 on 7 and 256 DF,  p-value: < 2.2e-16
plot_model(weight.curv, type = "pred", terms = c("day.num"),
           colors = kittens.coded$color)

Adding interaction terms for day number and the squared term with Cat estimates a curvilinear trajectory for each cat.

weight.curv.int<-lm(Weight ~ day.num*Cat + I(day.num^2)*Cat + Pre_post.eating +  Cat, kittens.coded)

summary(weight.curv.int)
## 
## Call:
## lm(formula = Weight ~ day.num * Cat + I(day.num^2) * Cat + Pre_post.eating + 
##     Cat, data = kittens.coded)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.075  -8.861   0.361   8.925  38.651 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               268.65596    6.80965  39.452  < 2e-16 ***
## day.num                    12.25007    0.88625  13.822  < 2e-16 ***
## CatGreen                   28.99428    9.48089   3.058  0.00247 ** 
## CatLavender                -4.22161    9.47675  -0.445  0.65637    
## CatMaroon                  -1.66805    9.47666  -0.176  0.86042    
## CatNo Tip                   2.27275    9.51064   0.239  0.81133    
## I(day.num^2)                0.14481    0.02483   5.833 1.69e-08 ***
## Pre_post.eatingPre-eating -22.57542    1.91595 -11.783  < 2e-16 ***
## day.num:CatGreen           -1.91283    1.25280  -1.527  0.12808    
## day.num:CatLavender        -1.34980    1.25251  -1.078  0.28223    
## day.num:CatMaroon           1.68849    1.25243   1.348  0.17883    
## day.num:CatNo Tip           3.21027    1.26743   2.533  0.01193 *  
## CatGreen:I(day.num^2)       0.16083    0.03510   4.582 7.30e-06 ***
## CatLavender:I(day.num^2)    0.21280    0.03510   6.063 4.94e-09 ***
## CatMaroon:I(day.num^2)     -0.01730    0.03509  -0.493  0.62246    
## CatNo Tip:I(day.num^2)     -0.03500    0.03617  -0.968  0.33412    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.37 on 248 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.994,  Adjusted R-squared:  0.9936 
## F-statistic:  2718 on 15 and 248 DF,  p-value: < 2.2e-16
plot_model(weight.curv.int, type = "pred", terms = c("day.num", "Cat"),
           colors = kittens.coded$color)