Reading in The Data
library(readxl)
cases <- read_xlsx("DVM3_Data.xlsx", "Data3")
control <- read_xlsx("DVM3_Data.xlsx","CaseControl2")
random <- read_xlsx("DVM3_Data.xlsx","Random")
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'require(foreign)## Loading required package: foreignlibrary(nnet)
require(ggplot2)## Loading required package: ggplot2require(reshape2)## Loading required package: reshape2library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(parameters)Descriptive Statistics for Case Population
glimpse(cases)## Rows: 162
## Columns: 19
## $ `Client Number` <chr> "1003644", "23322SU", "1019545", "30931SU", "10…
## $ `Patient Name` <chr> "Oakley", "Sasha", "Jackie", "Sky", "Cindy", "S…
## $ Patient_Breed <chr> "American Bulldog", "American Staffordshire Ter…
## $ AKC <chr> "Foundation Stock Service", "Terrier Group", "T…
## $ Breed_Category <chr> "Large", "Large", "Large", "Large", "Large", "L…
## $ `Post Code` <dbl> 2040, 2050, 2204, 2220, 2041, 2043, 2043, 2043,…
## $ Region <chr> "City", "City", "Inner West Sydney", "South Syd…
## $ Age_in_Months <dbl> 14, 79, 50, 27, 100, 32, 55, 66, 30, 12, 45, 20…
## $ Age_Category <chr> "<= 2 Y", "5-8 Y", "2-5 Y", "2-5 Y", ">= 8 Y", …
## $ Sex <chr> "M", "M", "F", "F", "F", "F", "F", "F", "M", "M…
## $ Neuter <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "N…
## $ Weight <dbl> 25.0, 29.5, 16.1, 23.4, 24.0, 21.6, 21.7, 19.3,…
## $ Weight_Category <chr> "> 20 kg", "> 20 kg", "10-20 kg", "> 20 kg", ">…
## $ `FB Identified` <chr> "Linear FB/Bone/Hair", "Rubber Toy", "Bone", "F…
## $ FB <chr> "Mixed", "Blunt", "Pointed", "Sharp", "Pointed"…
## $ Sock_Ingestion <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Samoyed <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N…
## $ Procedure <chr> "Surgery", "Surgery", "Surgery", "Surgery", "En…
## $ `Surgery vs Endoscopy` <chr> "Surgery - Gastrotomy + 2 Enterotomies (Duodenu…cases$`Sock_Ingestion`<- as.factor(cases$`Sock_Ingestion`)
cases$Samoyed <- as.factor(cases$Samoyed)
cases$`FB`<- as.factor(cases$`FB`)
cases$Sex <- as.factor(cases$Sex)
cases$Neuter <- as.factor(cases$Neuter)
cases$Age_Category <- as.factor(cases$Age_Category)
cases$Weight_Category <- as.factor(cases$Weight_Category)
cases$Breed_Category <- as.factor(cases$AKC)
cases$Procedure <- as.factor(cases$Procedure)
b <- table(cases$AKC)
b##
## Foundation Stock Service Herding Group Hound Group
## 3 4 13
## Mixed Breed Non-sporting Group Sporting Group
## 10 15 43
## Terrier Group Toy Group Working Group
## 40 21 13bpercent <- prop.table(b)*100
bpercent ##
## Foundation Stock Service Herding Group Hound Group
## 1.851852 2.469136 8.024691
## Mixed Breed Non-sporting Group Sporting Group
## 6.172840 9.259259 26.543210
## Terrier Group Toy Group Working Group
## 24.691358 12.962963 8.024691s <- table(cases$Sock_Ingestion)
s##
## 0 1
## 145 17spercent <- prop.table(s)*100
spercent ##
## 0 1
## 89.50617 10.49383samoyed <- table(cases$Samoyed)
samoyed ##
## N Y
## 157 5f <- table(cases$FB)
f##
## Blunt Linear Mixed Pointed Sharp
## 41 57 11 44 9fpercent <- prop.table(f)*100
fpercent ##
## Blunt Linear Mixed Pointed Sharp
## 25.308642 35.185185 6.790123 27.160494 5.555556g <- table(cases$Sex)
g##
## F M
## 61 101gpercent <- prop.table(g)*100
gpercent ##
## F M
## 37.65432 62.34568n <- table(cases$Neuter)
n ##
## N Y
## 30 132npercent <- prop.table(n)*100
npercent##
## N Y
## 18.51852 81.48148w <- table(cases$Weight_Category)
w ##
## < 5 kg > 20 kg 10-20 kg 5-10 kg
## 9 71 43 38wpercent <- prop.table(w)*100
wpercent##
## < 5 kg > 20 kg 10-20 kg 5-10 kg
## 5.590062 44.099379 26.708075 23.602484a <- table(cases$Age_Category)
a ##
## <= 2 Y >= 8 Y 2-5 Y 5-8 Y
## 56 38 46 22apercent <- prop.table(a)*100
apercent ##
## <= 2 Y >= 8 Y 2-5 Y 5-8 Y
## 34.56790 23.45679 28.39506 13.58025p <- table(cases$Procedure)
p##
## Endoscopy Surgery
## 54 108ppercent <- prop.table(p)*100
ppercent ##
## Endoscopy Surgery
## 33.33333 66.66667bs <- table(cases$Sock_Ingestion,cases$Breed_Category)
table(cases$Sock_Ingestion,cases$Weight_Category)##
## < 5 kg > 20 kg 10-20 kg 5-10 kg
## 0 8 62 39 36
## 1 1 9 4 2table(cases$Sock_Ingestion,cases$Sex)##
## F M
## 0 55 90
## 1 6 11rb <- table(random$AKC)
rb##
## Foundation Stock Service Herding Group Hound Group
## 13 29 40
## Mixed Breed Non-sporting Group Sporting Group
## 36 30 41
## Terrier Group Toy Group Working Group
## 63 56 21rbpercent <- prop.table(rb)*100
rbpercent##
## Foundation Stock Service Herding Group Hound Group
## 3.951368 8.814590 12.158055
## Mixed Breed Non-sporting Group Sporting Group
## 10.942249 9.118541 12.462006
## Terrier Group Toy Group Working Group
## 19.148936 17.021277 6.382979median(random$Age_in_Months)## [1] 64ra <- table(random$Age_Category)
rapercent <- prop.table(ra)*100
rapercent##
## < = 2 Y > = 8 Y 2-5 Y 5-8 Y
## 29.48328 31.30699 17.32523 21.88450median(random$Weight)## [1] 14.5rw <- table(random$Weight_Category)
rwpercent <- prop.table(rw)*100
rwpercent ##
## <5 kg >20 kg 10-20 kg 5-10 kg
## 11.85410 40.12158 23.10030 24.92401rs <- table(random$Sex)
rs##
## F M
## 133 196rspercent <- prop.table(rs)*100
rspercent##
## F M
## 40.42553 59.57447rn <- table(random$Neuter)
rn##
## N Y
## 111 218rnpercent <- prop.table(rn)*100
rnpercent##
## N Y
## 33.7386 66.2614Descriptive Statistics for Case + Random (Case-Control) Population Group
glimpse(control)## Rows: 491
## Columns: 14
## $ `Client Number` <chr> "31604SU", "1003272", "23808SU", "1003208", "1000066",…
## $ `Patient Name` <chr> "Coco", "Tiger", "Mr Bean", "Cash", "Sarge", "Newton",…
## $ `Patient Breed` <chr> "Alaskan Malamute", "American Bulldog", "American Staf…
## $ AKC <chr> "Working Group", "Foundation Stock Service", "Terrier …
## $ `Sex + Neuter` <chr> "F", "M", "M", "M", "M", "M", "F", "M", "M", "M", "M",…
## $ Sex <chr> "F", "M", "M", "M", "M", "M", "F", "M", "M", "M", "M",…
## $ Neuter <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ Weight <dbl> 39.60, 33.30, 30.20, 32.60, 21.00, 29.50, 26.40, 24.50…
## $ Weight_Category <chr> "> 20 kg", "> 20 kg", "> 20 kg", "> 20 kg", "> 20 kg",…
## $ Age_in_Months <dbl> 88, 90, 19, 132, 49, 11, 72, 72, 96, 2, 16, 160, 2, 36…
## $ Age_Category <chr> "5-8 Y", "5-8 Y", "< = 2 Y", "> = 8 Y", "2-5 Y", "< = …
## $ `FB (Y or N)` <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ FB <chr> "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0",…
## $ Samoyed <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …control$Breed <- as.factor(control$AKC)
control$Sex <- as.factor(control$Sex)
control$Neuter <- as.factor(control$Neuter)
control$Weight_Category <- as.factor(control$Weight_Category)
control$Age_Category <- as.factor(control$Age_Category)
control$FB <- as.factor(control$FB)
control$Samoyed <- as.factor(control$Samoyed)
sex <- table(control$Sex)
sex##
## F M
## 194 297sexpercent <- prop.table(sex)*100
sexpercent##
## F M
## 39.5112 60.4888neuter <- table(control$Neuter)
neuter##
## N Y
## 141 350neuterpercent <- prop.table(neuter)*100
neuterpercent##
## N Y
## 28.7169 71.2831foreign <- table(control$FB)
foreign ##
## 0 1
## 329 162foreignpercent <- prop.table(foreign)*100
foreignpercent ##
## 0 1
## 67.00611 32.99389breed <- table(control$AKC)
breed ##
## Foundation Stock Service Herding Group Hound Group
## 13 36 53
## Mixed Breed Non-sporting Group Sporting Group
## 47 44 83
## Terrier Group Toy Group Working Group
## 103 77 35breedpercent <- prop.table(breed)*100
breedpercent ##
## Foundation Stock Service Herding Group Hound Group
## 2.647658 7.331976 10.794297
## Mixed Breed Non-sporting Group Sporting Group
## 9.572301 8.961303 16.904277
## Terrier Group Toy Group Working Group
## 20.977597 15.682281 7.128310weight <- table(control$Weight_Category)
weight ##
## < 5 kg > 20 kg 10-20 kg 5-10 kg
## 48 204 119 120weightpercent <- prop.table(weight)*100
weightpercent ##
## < 5 kg > 20 kg 10-20 kg 5-10 kg
## 9.775967 41.547862 24.236253 24.439919age <- table(control$Age_Category)
age ##
## < = 2 Y > = 8 Y 2-5 Y 5-8 Y
## 153 141 103 94agepercent <- prop.table(age)*100
agepercent ##
## < = 2 Y > = 8 Y 2-5 Y 5-8 Y
## 31.1609 28.7169 20.9776 19.1446Univariable Logistic Regression (Case-Control) Sex
SexControlModel <- glm(FB ~ Sex, data = control, family = binomial(link = "logit"))
summary(SexControlModel)##
## Call:
## glm(formula = FB ~ Sex, family = binomial(link = "logit"), data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7795 0.1546 -5.041 4.64e-07 ***
## SexM 0.1165 0.1973 0.590 0.555
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 622.37 on 489 degrees of freedom
## AIC: 626.37
##
## Number of Fisher Scoring iterations: 4anova(SexControlModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## Sex 1 0.34984 489 622.37 0.5542odds_ratios_s <- exp(coef(SexControlModel))
print(odds_ratios_s)## (Intercept) SexM
## 0.4586466 1.1235363conf_int_s <- confint(SexControlModel)## Waiting for profiling to be done...conf_int_s_odds_ratio <- exp(conf_int_s)
print(conf_int_s_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.3365974 0.6178629
## SexM 0.7646804 1.6584872Neuter
NeuterControlModel <- glm (FB ~ Neuter, data = control, family = binomial(link = "logit"))
summary(NeuterControlModel) ##
## Call:
## glm(formula = FB ~ Neuter, family = binomial(link = "logit"),
## data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3083 0.2058 -6.358 2.04e-10 ***
## NeuterY 0.8066 0.2335 3.455 0.00055 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 609.82 on 489 degrees of freedom
## AIC: 613.82
##
## Number of Fisher Scoring iterations: 4anova(NeuterControlModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## Neuter 1 12.905 489 609.82 0.0003278 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios_n <- exp(coef(NeuterControlModel))
print(odds_ratios_n)## (Intercept) NeuterY
## 0.2702703 2.2403670conf_int_n <- confint(NeuterControlModel)## Waiting for profiling to be done...conf_int_n_odds_ratio <- exp(conf_int_n)
print(conf_int_n_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.1774972 0.3988573
## NeuterY 1.4326442 3.5868737Breed
control$FB <- as.factor(control$FB)
control$AKC <- as.factor(control$AKC)
p.adjust(control$AKC, method = "bonferroni")## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [482] 1 1 1 1 1 1 1 1 1 1control$BreedControl <- relevel(control$AKC, ref = "Mixed Breed")
BreedControlModel <- glm(FB ~ AKC, data = control, family = binomial(link = "logit"))
summary(BreedControlModel)##
## Call:
## glm(formula = FB ~ AKC, family = binomial(link = "logit"), data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.20397 0.65828 -1.829 0.0674 .
## AKCHerding Group -0.87547 0.84532 -1.036 0.3004
## AKCHound Group 0.08004 0.73161 0.109 0.9129
## AKCMixed Breed 0.01835 0.74298 0.025 0.9803
## AKCNon-sporting Group 0.44183 0.73355 0.602 0.5470
## AKCSporting Group 1.22807 0.69393 1.770 0.0768 .
## AKCTerrier Group 0.74972 0.68863 1.089 0.2763
## AKCToy Group 0.22314 0.70626 0.316 0.7520
## AKCWorking Group 0.79851 0.74322 1.074 0.2827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 594.41 on 482 degrees of freedom
## AIC: 612.41
##
## Number of Fisher Scoring iterations: 4anova(BreedControlModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## AKC 8 28.311 482 594.41 0.0004186 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios_b <- exp(coef(BreedControlModel))
print(odds_ratios_b)## (Intercept) AKCHerding Group AKCHound Group
## 0.3000000 0.4166667 1.0833333
## AKCMixed Breed AKCNon-sporting Group AKCSporting Group
## 1.0185185 1.5555556 3.4146341
## AKCTerrier Group AKCToy Group AKCWorking Group
## 2.1164021 1.2500000 2.2222222conf_int_b <- confint(BreedControlModel)## Waiting for profiling to be done...conf_int_b_odds_ratio <- exp(conf_int_b)
print(conf_int_b_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.06725892 0.9806492
## AKCHerding Group 0.07822490 2.4060611
## AKCHound Group 0.27915838 5.3690221
## AKCMixed Breed 0.25517717 5.1303105
## AKCNon-sporting Group 0.40000047 7.7394427
## AKCSporting Group 0.96443288 16.0438608
## AKCTerrier Group 0.60442615 9.8634404
## AKCToy Group 0.34236900 5.9733283
## AKCWorking Group 0.56050013 11.2332467Weight
control$Weight_Category <- relevel(control$Weight_Category, ref = "10-20 kg")
WeightControlModel <- glm(FB ~ Weight_Category, data = control, family = binomial(link="logit"))
summary(WeightControlModel)##
## Call:
## glm(formula = FB ~ Weight_Category, family = binomial(link = "logit"),
## data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5695 0.1908 -2.985 0.00284 **
## Weight_Category< 5 kg -0.8968 0.4161 -2.155 0.03115 *
## Weight_Category> 20 kg -0.0366 0.2406 -0.152 0.87907
## Weight_Category5-10 kg -0.1996 0.2737 -0.729 0.46588
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 616.76 on 487 degrees of freedom
## AIC: 624.76
##
## Number of Fisher Scoring iterations: 4anova(WeightControlModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## Weight_Category 3 5.9639 487 616.76 0.1134odds_ratios_w <- exp(coef(WeightControlModel))
print(odds_ratios_w)## (Intercept) Weight_Category< 5 kg Weight_Category> 20 kg
## 0.5657895 0.4078712 0.9640592
## Weight_Category5-10 kg
## 0.8190584conf_int_w <- confint(WeightControlModel)## Waiting for profiling to be done...conf_int_w_odds_ratio <- exp(conf_int_w)
print(conf_int_w_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.3862924 0.8179947
## Weight_Category< 5 kg 0.1715035 0.8905444
## Weight_Category> 20 kg 0.6027101 1.5500830
## Weight_Category5-10 kg 0.4777488 1.3999556Age
control$Age_Category <- relevel(control$Age_Category, ref = "5-8 Y")
AgeControlModel <- glm(FB ~ Age_Category, data = control, family = binomial(link = "logit"))
summary(AgeControlModel)##
## Call:
## glm(formula = FB ~ Age_Category, family = binomial(link = "logit"),
## data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1856 0.2436 -4.867 1.13e-06 ***
## Age_Category< = 2 Y 0.6363 0.2958 2.151 0.03149 *
## Age_Category> = 8 Y 0.1885 0.3088 0.610 0.54164
## Age_Category2-5 Y 0.9712 0.3140 3.093 0.00198 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 609.23 on 487 degrees of freedom
## AIC: 617.23
##
## Number of Fisher Scoring iterations: 4anova(AgeControlModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## Age_Category 3 13.496 487 609.23 0.003678 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios <- exp(coef(AgeControlModel))
print(odds_ratios)## (Intercept) Age_Category< = 2 Y Age_Category> = 8 Y Age_Category2-5 Y
## 0.3055556 1.8894096 1.2074139 2.6411483conf_int_a <- confint(AgeControlModel)## Waiting for profiling to be done...conf_int_a_odds_ratio <- exp(conf_int_a)
print(conf_int_a_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.1852797 0.4838347
## Age_Category< = 2 Y 1.0690256 3.4224872
## Age_Category> = 8 Y 0.6634174 2.2359744
## Age_Category2-5 Y 1.4405155 4.9521850SamoyedControlModel <- glm(FB ~ Samoyed, data = control, family = binomial(link = "logit"))
summary(SamoyedControlModel)##
## Call:
## glm(formula = FB ~ Samoyed, family = binomial(link = "logit"),
## data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.76913 0.09812 -7.839 4.56e-15 ***
## Samoyed1 3.07172 1.05317 2.917 0.00354 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 606.06 on 489 degrees of freedom
## AIC: 610.06
##
## Number of Fisher Scoring iterations: 4anova(SamoyedControlModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## Samoyed 1 16.661 489 606.06 4.469e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios_s <- exp(coef(SamoyedControlModel))
print(odds_ratios_s)## (Intercept) Samoyed1
## 0.4634146 21.5789452conf_int_s <- confint(SamoyedControlModel)## Waiting for profiling to be done...conf_int_s_odds_ratio <- exp(conf_int_s)
print(conf_int_s_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.3814148 0.5604979
## Samoyed1 4.0793155 397.8188587The variables that have a significant impact on FB ingestion include Age Category, Weight Category and Neuter status.
Multivariable Logistic Regression (Case-Control)
null_model <- glm(FB ~ 1, data = control, family = binomial)
full_model <- glm(FB ~ Age_Category + Neuter + AKC, data = control, family = binomial(link = "logit"))
stepwise_model <- step(null_model, scope = list(lower = null_model, upper = full_model), direction = "forward")## Start: AIC=624.72
## FB ~ 1
##
## Df Deviance AIC
## + AKC 8 594.41 612.41
## + Neuter 1 609.82 613.82
## + Age_Category 3 609.23 617.23
## <none> 622.72 624.72
##
## Step: AIC=612.41
## FB ~ AKC
##
## Df Deviance AIC
## + Neuter 1 579.58 599.58
## + Age_Category 3 580.16 604.16
## <none> 594.41 612.41
##
## Step: AIC=599.58
## FB ~ AKC + Neuter
##
## Df Deviance AIC
## + Age_Category 3 558.95 584.95
## <none> 579.58 599.58
##
## Step: AIC=584.95
## FB ~ AKC + Neuter + Age_Categoryanova(null_model, full_model, test = "LRT")## Analysis of Deviance Table
##
## Model 1: FB ~ 1
## Model 2: FB ~ Age_Category + Neuter + AKC
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 490 622.72
## 2 478 558.95 12 63.773 4.592e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(full_model)##
## Call:
## glm(formula = FB ~ Age_Category + Neuter + AKC, family = binomial(link = "logit"),
## data = control)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.86083 0.78248 -3.656 0.000256 ***
## Age_Category< = 2 Y 1.05147 0.33079 3.179 0.001479 **
## Age_Category> = 8 Y 0.17018 0.32395 0.525 0.599366
## Age_Category2-5 Y 1.11517 0.33566 3.322 0.000893 ***
## NeuterY 1.15040 0.26286 4.376 1.21e-05 ***
## AKCHerding Group -0.82744 0.86710 -0.954 0.339948
## AKCHound Group 0.25207 0.75786 0.333 0.739434
## AKCMixed Breed -0.08224 0.76747 -0.107 0.914663
## AKCNon-sporting Group 0.53127 0.75839 0.701 0.483595
## AKCSporting Group 1.35423 0.71782 1.887 0.059214 .
## AKCTerrier Group 1.08417 0.71363 1.519 0.128702
## AKCToy Group 0.46313 0.72940 0.635 0.525465
## AKCWorking Group 0.95347 0.76824 1.241 0.214562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 622.72 on 490 degrees of freedom
## Residual deviance: 558.95 on 478 degrees of freedom
## AIC: 584.95
##
## Number of Fisher Scoring iterations: 4anova(full_model, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: FB
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 490 622.72
## Age_Category 3 13.496 487 609.23 0.003678 **
## Neuter 1 18.329 486 590.90 1.859e-05 ***
## AKC 8 31.949 478 558.95 9.513e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios1 <- exp(coef(full_model))
print(odds_ratios1)## (Intercept) Age_Category< = 2 Y Age_Category> = 8 Y
## 0.05722123 2.86185902 1.18551541
## Age_Category2-5 Y NeuterY AKCHerding Group
## 3.05008255 3.15945513 0.43716622
## AKCHound Group AKCMixed Breed AKCNon-sporting Group
## 1.28668110 0.92104929 1.70109735
## AKCSporting Group AKCTerrier Group AKCToy Group
## 3.87378845 2.95698860 1.58904030
## AKCWorking Group
## 2.59470244conf_int <- confint(full_model)## Waiting for profiling to be done...conf_int_odds_ratio <- exp(conf_int)
print(conf_int_odds_ratio)## 2.5 % 97.5 %
## (Intercept) 0.01054579 0.2424444
## Age_Category< = 2 Y 1.51449071 5.5577506
## Age_Category> = 8 Y 0.63188195 2.2593605
## Age_Category2-5 Y 1.59567401 5.9697455
## NeuterY 1.91093474 5.3676525
## AKCHerding Group 0.07868058 2.6176252
## AKCHound Group 0.31297135 6.6388155
## AKCMixed Breed 0.21836027 4.8137352
## AKCNon-sporting Group 0.41397736 8.7915891
## AKCSporting Group 1.03671199 18.8590227
## AKCTerrier Group 0.79985817 14.3133774
## AKCToy Group 0.41366307 7.8665661
## AKCWorking Group 0.61969597 13.6375337Age Category For animals in the age category “8 years or older” compared to the baseline age category, the log odds of ‘FB’ decrease by 0.7892. In terms of odds ratio, the odds of ‘FB’ decrease by 0.453 in animals aged 8 yeras or older compared to the baseline age category.
Proportion Tables for Sock Ingestion in Case Group Breed
table(cases$Breed_Category, cases$Sock_Ingestion)##
## 0 1
## Foundation Stock Service 3 0
## Herding Group 4 0
## Hound Group 13 0
## Mixed Breed 8 2
## Non-sporting Group 14 1
## Sporting Group 36 7
## Terrier Group 40 0
## Toy Group 20 1
## Working Group 7 6Poodle Mix - Sporting Group
prop.test(c(3,6),c(12,37))## Warning in prop.test(c(3, 6), c(12, 37)): Chi-squared approximation may be
## incorrect##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(3, 6) out of c(12, 37)
## X-squared = 0.064454, df = 1, p-value = 0.7996
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.2396084 0.4152841
## sample estimates:
## prop 1 prop 2
## 0.2500000 0.1621622Poodle Mix - Toy Group
prop.test (c(3,2),c(12,22))## Warning in prop.test(c(3, 2), c(12, 22)): Chi-squared approximation may be
## incorrect##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(3, 2) out of c(12, 22)
## X-squared = 0.55512, df = 1, p-value = 0.4562
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.1781648 0.4963466
## sample estimates:
## prop 1 prop 2
## 0.25000000 0.09090909Poodle Mix - Working Group
prop.test(c(3,4),c(12,10))## Warning in prop.test(c(3, 4), c(12, 10)): Chi-squared approximation may be
## incorrect##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(3, 4) out of c(12, 10)
## X-squared = 0.085556, df = 1, p-value = 0.7699
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.6318177 0.3318177
## sample estimates:
## prop 1 prop 2
## 0.25 0.40Sporting Group- Toy Group
prop.test(c(6,2),c(37,22))## Warning in prop.test(c(6, 2), c(37, 22)): Chi-squared approximation may be
## incorrect##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(6, 2) out of c(37, 22)
## X-squared = 0.1443, df = 1, p-value = 0.704
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.1339159 0.2764220
## sample estimates:
## prop 1 prop 2
## 0.16216216 0.09090909Sporting Group - Working Group
prop.test(c(6,4),c(37,10))## Warning in prop.test(c(6, 4), c(37, 10)): Chi-squared approximation may be
## incorrect##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(6, 4) out of c(37, 10)
## X-squared = 1.4283, df = 1, p-value = 0.232
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.6273897 0.1517140
## sample estimates:
## prop 1 prop 2
## 0.1621622 0.4000000There is no significant difference in the proportion of animals that have ingested sock and not ingested sock between any of the following Breed Groups - Poodle Mix, Toy Group, Working Group, Sporting Group.
SamoyedCaseModel <- glm(Sock_Ingestion ~ Samoyed, data = cases, family = binomial(link = "logit"))
summary(SamoyedCaseModel)##
## Call:
## glm(formula = Sock_Ingestion ~ Samoyed, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3238 0.2800 -8.298 < 2e-16 ***
## SamoyedY 2.7293 0.9549 2.858 0.00426 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.80 on 161 degrees of freedom
## Residual deviance: 101.12 on 160 degrees of freedom
## AIC: 105.12
##
## Number of Fisher Scoring iterations: 5anova(SamoyedCaseModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 161 108.80
## Samoyed 1 7.675 160 101.12 0.005599 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios2 <- exp(coef(SamoyedCaseModel))
print(odds_ratios2)## (Intercept) SamoyedY
## 0.0979021 15.3214286conf_int2 <- confint(SamoyedCaseModel)## Waiting for profiling to be done...conf_int_odds_ratio2 <- exp(conf_int2)
print(conf_int_odds_ratio2)## 2.5 % 97.5 %
## (Intercept) 0.0540154 0.1632807
## SamoyedY 2.3565726 124.0634703AgeCaseModel <- glm(Sock_Ingestion ~ Age_Category, data = cases, family = binomial(link = "logit"))
anova(AgeCaseModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 161 108.80
## Age_Category 3 6.6696 158 102.13 0.08321 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(AgeCaseModel)##
## Call:
## glm(formula = Sock_Ingestion ~ Age_Category, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6529 0.3639 -4.543 5.55e-06 ***
## Age_Category>= 8 Y -0.4871 0.6417 -0.759 0.448
## Age_Category2-5 Y -0.6985 0.6373 -1.096 0.273
## Age_Category5-8 Y -16.9131 1390.6313 -0.012 0.990
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.80 on 161 degrees of freedom
## Residual deviance: 102.13 on 158 degrees of freedom
## AIC: 110.13
##
## Number of Fisher Scoring iterations: 17odds_ratio_age <- exp(coef(AgeCaseModel))
print(odds_ratio_age)## (Intercept) Age_Category>= 8 Y Age_Category2-5 Y Age_Category5-8 Y
## 1.914894e-01 6.143791e-01 4.973545e-01 4.515587e-08conf_int_age <- confint(AgeCaseModel)## Waiting for profiling to be done...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredconf_int_age2 <- exp(conf_int_age)
print(conf_int_age2)## 2.5 % 97.5 %
## (Intercept) 0.08776388 3.718334e-01
## Age_Category>= 8 Y 0.15604540 2.056579e+00
## Age_Category2-5 Y 0.12718147 1.648788e+00
## Age_Category5-8 Y NA 2.902321e+31SexCaseModel <- glm(Sock_Ingestion ~ Sex, data = cases, family = binomial(link = "logit"))
summary(SexCaseModel)##
## Call:
## glm(formula = Sock_Ingestion ~ Sex, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2156 0.4299 -5.154 2.55e-07 ***
## SexM 0.1137 0.5356 0.212 0.832
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.80 on 161 degrees of freedom
## Residual deviance: 108.75 on 160 degrees of freedom
## AIC: 112.75
##
## Number of Fisher Scoring iterations: 4anova(SexCaseModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 161 108.80
## Sex 1 0.045418 160 108.75 0.8312odds_ratio_sex <- exp(coef(SexCaseModel))
print(odds_ratio_sex)## (Intercept) SexM
## 0.1090909 1.1203703conf_int_sex <- confint(SexCaseModel)## Waiting for profiling to be done...conf_int_sex2 <- exp(conf_int_sex)
print(conf_int_sex2)## 2.5 % 97.5 %
## (Intercept) 0.04199906 0.2335333
## SexM 0.40241824 3.4100146NeuterCaseModel <- glm(Sock_Ingestion ~ Neuter, data = cases, family = binomial(link = "logit"))
summary(NeuterCaseModel)##
## Call:
## glm(formula = Sock_Ingestion ~ Neuter, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3863 0.4564 -3.037 0.00239 **
## NeuterY -1.0116 0.5545 -1.824 0.06812 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.80 on 161 degrees of freedom
## Residual deviance: 105.75 on 160 degrees of freedom
## AIC: 109.75
##
## Number of Fisher Scoring iterations: 5anova(NeuterCaseModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 161 108.80
## Neuter 1 3.0503 160 105.75 0.08072 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratio_neuter <- exp(coef(NeuterCaseModel))
print(odds_ratio_neuter)## (Intercept) NeuterY
## 0.2500000 0.3636364conf_int_neuter <- confint(NeuterCaseModel)## Waiting for profiling to be done...conf_int_neuter2 <- exp(conf_int_neuter)
print(conf_int_neuter2)## 2.5 % 97.5 %
## (Intercept) 0.09261171 0.5731953
## NeuterY 0.12525514 1.1409732BreedCaseModel <- glm(Sock_Ingestion ~ Breed_Category, data = cases, family = binomial(link = "logit"))
summary(BreedCaseModel)##
## Call:
## glm(formula = Sock_Ingestion ~ Breed_Category, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.957e+01 6.209e+03 -0.003 0.997
## Breed_CategoryHerding Group -2.791e-08 8.214e+03 0.000 1.000
## Breed_CategoryHound Group -2.791e-08 6.888e+03 0.000 1.000
## Breed_CategoryMixed Breed 1.818e+01 6.209e+03 0.003 0.998
## Breed_CategoryNon-sporting Group 1.693e+01 6.209e+03 0.003 0.998
## Breed_CategorySporting Group 1.793e+01 6.209e+03 0.003 0.998
## Breed_CategoryTerrier Group -2.792e-08 6.437e+03 0.000 1.000
## Breed_CategoryToy Group 1.657e+01 6.209e+03 0.003 0.998
## Breed_CategoryWorking Group 1.941e+01 6.209e+03 0.003 0.998
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.799 on 161 degrees of freedom
## Residual deviance: 81.549 on 153 degrees of freedom
## AIC: 99.549
##
## Number of Fisher Scoring iterations: 18anova(BreedCaseModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 161 108.799
## Breed_Category 8 27.251 153 81.549 0.0006399 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratio_breed <- exp(coef(BreedCaseModel))
print(odds_ratio_breed)## (Intercept) Breed_CategoryHerding Group
## 3.181006e-09 1.000000e+00
## Breed_CategoryHound Group Breed_CategoryMixed Breed
## 1.000000e+00 7.859150e+07
## Breed_CategoryNon-sporting Group Breed_CategorySporting Group
## 2.245471e+07 6.112672e+07
## Breed_CategoryTerrier Group Breed_CategoryToy Group
## 1.000000e+00 1.571830e+07
## Breed_CategoryWorking Group
## 2.694566e+08conf_int_breed <- confint(BreedCaseModel)## Waiting for profiling to be done...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredconf_int_breed2 <- exp(conf_int_breed)
print(conf_int_breed2)## 2.5 % 97.5 %
## (Intercept) NA Inf
## Breed_CategoryHerding Group 5.496609e-104 1.819270e+103
## Breed_CategoryHound Group 1.969412e-73 6.343726e+65
## Breed_CategoryMixed Breed 7.737750e-165 NA
## Breed_CategoryNon-sporting Group 2.210780e-165 NA
## Breed_CategorySporting Group 0.000000e+00 NA
## Breed_CategoryTerrier Group 0.000000e+00 2.199771e+165
## Breed_CategoryToy Group 1.547546e-165 NA
## Breed_CategoryWorking Group 0.000000e+00 NAWeightCaseModel <- glm(Sock_Ingestion ~ Weight_Category, data = cases, family = binomial(link = "logit"))
summary(WeightCaseModel)##
## Call:
## glm(formula = Sock_Ingestion ~ Weight_Category, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0794 1.0607 -1.961 0.0499 *
## Weight_Category> 20 kg 0.1495 1.1190 0.134 0.8937
## Weight_Category10-20 kg -0.1978 1.1835 -0.167 0.8672
## Weight_Category5-10 kg -0.8109 1.2856 -0.631 0.5282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 104.24 on 160 degrees of freedom
## Residual deviance: 102.55 on 157 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 110.55
##
## Number of Fisher Scoring iterations: 5anova(WeightCaseModel, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 160 104.24
## Weight_Category 3 1.6861 157 102.55 0.64odds_ratio_weight <- exp(coef(WeightCaseModel))
print(odds_ratio_weight)## (Intercept) Weight_Category> 20 kg Weight_Category10-20 kg
## 0.1250000 1.1612903 0.8205128
## Weight_Category5-10 kg
## 0.4444444conf_int_weight <- confint(WeightCaseModel)## Waiting for profiling to be done...conf_int_weight2 <- exp(conf_int_weight)
print(conf_int_weight2)## 2.5 % 97.5 %
## (Intercept) 0.006737736 0.6811608
## Weight_Category> 20 kg 0.179615970 22.8403716
## Weight_Category10-20 kg 0.102849260 17.1750199
## Weight_Category5-10 kg 0.037755807 10.2469917Model2 <- glm(Sock_Ingestion ~ Samoyed + Neuter, data = cases, family = binomial(link = "logit"))
summary(Model2)##
## Call:
## glm(formula = Sock_Ingestion ~ Samoyed + Neuter, family = binomial(link = "logit"),
## data = cases)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3863 0.4564 -3.037 0.00239 **
## SamoyedY 3.1051 0.9832 3.158 0.00159 **
## NeuterY -1.3134 0.5846 -2.247 0.02466 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.799 on 161 degrees of freedom
## Residual deviance: 96.475 on 159 degrees of freedom
## AIC: 102.48
##
## Number of Fisher Scoring iterations: 5anova(Model2, test = 'Chisq')## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Sock_Ingestion
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 161 108.799
## Samoyed 1 7.6750 160 101.124 0.005599 **
## Neuter 1 4.6488 159 96.475 0.031075 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1odds_ratios3 <- exp(coef(Model2))
print(odds_ratios3)## (Intercept) SamoyedY NeuterY
## 0.2500000 22.3125000 0.2689076conf_int <- confint(Model2)## Waiting for profiling to be done...conf_int_odds_ratio2 <- exp(conf_int)
print(conf_int_odds_ratio2)## 2.5 % 97.5 %
## (Intercept) 0.09261171 0.5731953
## SamoyedY 3.27312320 189.5335373
## NeuterY 0.08552744 0.8814555Nominal Logistic Regression
NomModel <- multinom(FB ~ Age_Category, data = cases)## # weights: 25 (16 variable)
## initial value 260.728942
## iter 10 value 217.940294
## iter 20 value 217.546146
## iter 30 value 217.538501
## final value 217.538488
## convergedsummary_model <- summary(NomModel)
print(summary_model)## Call:
## multinom(formula = FB ~ Age_Category, data = cases)
##
## Coefficients:
## (Intercept) Age_Category>= 8 Y Age_Category2-5 Y Age_Category5-8 Y
## Linear 0.3022610 0.5732048 -0.3712425 0.3909072
## Mixed -0.8873324 0.3765477 -1.8204483 -13.5997620
## Pointed -0.8872864 2.0504602 0.7441958 1.5804403
## Sharp -2.1401420 1.2239588 0.5307039 1.4470110
##
## Std. Errors:
## (Intercept) Age_Category>= 8 Y Age_Category2-5 Y Age_Category5-8 Y
## Linear 0.3198459 0.6209992 0.4903035 0.6908713
## Mixed 0.4490906 0.8573281 1.1260962 699.5283419
## Pointed 0.4490832 0.6813068 0.5875938 0.7593940
## Sharp 0.7475657 1.1219684 0.9792137 1.1440521
##
## Residual Deviance: 435.077
## AIC: 467.077coefficients <- coef(NomModel)
std_errors <- summary_model$standard.errors
z_values <- coefficients/std_errors
p_values <- 2*(1 - pnorm(abs(z_values)))
print(p_values)## (Intercept) Age_Category>= 8 Y Age_Category2-5 Y Age_Category5-8 Y
## Linear 0.344648233 0.355988293 0.4489488 0.57151774
## Mixed 0.048172572 0.660508864 0.1059644 0.98448904
## Pointed 0.048180540 0.002615929 0.2053291 0.03741688
## Sharp 0.004198963 0.275315547 0.5878395 0.20593868odds_ratio <- exp(summary(NomModel)$coefficients)
data.frame(t(odds_ratio))## Linear Mixed Pointed Sharp
## (Intercept) 1.3529143 4.117527e-01 0.4117716 0.1176381
## Age_Category>= 8 Y 1.7739430 1.457245e+00 7.7714769 3.4006234
## Age_Category2-5 Y 0.6898767 1.619531e-01 2.1047481 1.7001286
## Age_Category5-8 Y 1.4783213 1.240790e-06 4.8570937 4.2503909‘Linear’ The odds of being classified as ‘Linear’ (versus the Blunt) increases by a factor of 1.4733 for individuals in the ‘<= 2 Y’ age category.
NomModel1 <- multinom(FB ~ Sex, data = cases)## # weights: 15 (8 variable)
## initial value 260.728942
## iter 10 value 225.281632
## final value 224.888646
## convergedsummary_model1 <- summary(NomModel1)
print(summary_model1)## Call:
## multinom(formula = FB ~ Sex, data = cases)
##
## Coefficients:
## (Intercept) SexM
## Linear 0.6061162 -0.41803147
## Mixed -1.3863626 0.09848846
## Pointed 0.3482686 -0.41970695
## Sharp -0.5390658 -2.13508508
##
## Std. Errors:
## (Intercept) SexM
## Linear 0.3588663 0.4379946
## Mixed 0.6455055 0.7590553
## Pointed 0.3770345 0.4622500
## Sharp 0.4755984 0.8721730
##
## Residual Deviance: 449.7773
## AIC: 465.7773coefficients1 <- coef(NomModel1)
std_errors1 <- summary_model1$standard.errors
z_values1 <- coefficients1/std_errors1
p_values1 <- 2*(1 - pnorm(abs(z_values1)))
print(p_values1)## (Intercept) SexM
## Linear 0.09122426 0.33987030
## Mixed 0.03173631 0.89676314
## Pointed 0.35563992 0.36389660
## Sharp 0.25702635 0.01436491odds_ratio1 <- exp(summary(NomModel1)$coefficients)
data.frame(t(odds_ratio1))## Linear Mixed Pointed Sharp
## (Intercept) 1.8332974 0.2499829 1.4166127 0.5832929
## SexM 0.6583415 1.1035017 0.6572394 0.1182345NomModel2 <- multinom(FB ~ Neuter, data = cases)## # weights: 15 (8 variable)
## initial value 260.728942
## iter 10 value 224.290623
## final value 223.997880
## convergedsummary_model2 <- summary(NomModel2)
print(summary_model2)## Call:
## multinom(formula = FB ~ Neuter, data = cases)
##
## Coefficients:
## (Intercept) NeuterY
## Linear 0.3675201 -0.04931483
## Mixed -0.8117062 -0.70816138
## Pointed -0.8111388 1.03434796
## Sharp -8.8552849 7.58672627
##
## Std. Errors:
## (Intercept) NeuterY
## Linear 0.4336088 0.4919307
## Mixed 0.6010334 0.7316710
## Pointed 0.6009154 0.6460215
## Sharp 27.9107770 27.9133271
##
## Residual Deviance: 447.9958
## AIC: 463.9958coefficients2 <- coef(NomModel2)
std_errors2 <- summary_model2$standard.errors
z_values2 <- coefficients2/std_errors2
p_values2 <- 2*(1 - pnorm(abs(z_values2)))
print(p_values2)## (Intercept) NeuterY
## Linear 0.3966695 0.9201478
## Mixed 0.1768500 0.3331100
## Pointed 0.1770677 0.1093538
## Sharp 0.7510379 0.7857790odds_ratio2 <- exp(summary(NomModel2)$coefficients)
data.frame(t(odds_ratio2))## Linear Mixed Pointed Sharp
## (Intercept) 1.4441488 0.4440997 0.4443517 1.426260e-04
## NeuterY 0.9518814 0.4925490 2.8132713 1.971848e+03NomModel3 <- multinom(FB ~ Breed_Category, data = cases)## # weights: 50 (36 variable)
## initial value 260.728942
## iter 10 value 207.110423
## iter 20 value 201.229980
## iter 30 value 200.314790
## iter 40 value 200.293775
## final value 200.293709
## convergedsummary_model3 <- summary(NomModel3)## Warning in sqrt(diag(vc)): NaNs producedprint(summary_model3)## Call:
## multinom(formula = FB ~ Breed_Category, data = cases)
##
## Coefficients:
## (Intercept) Breed_CategoryHerding Group Breed_CategoryHound Group
## Linear -21.38613 22.07911 22.89044
## Mixed 21.13865 -45.96189 -21.83138
## Pointed -19.86609 -38.16923 19.17336
## Sharp 20.44541 -20.44556 -39.36754
## Breed_CategoryMixed Breed Breed_CategoryNon-sporting Group
## Linear 21.89693 21.38613
## Mixed -61.58694 -22.23736
## Pointed 19.46057 20.71335
## Sharp -45.61159 -21.54400
## Breed_CategorySporting Group Breed_CategoryTerrier Group
## Linear 21.13482 21.60926
## Mixed -22.93042 -22.52496
## Pointed 19.05516 20.49468
## Sharp -46.20947 -20.91543
## Breed_CategoryToy Group Breed_CategoryWorking Group
## Linear 21.72261 23.33205
## Mixed -60.63515 -20.44548
## Pointed 20.33610 20.96472
## Sharp -22.05480 -33.93285
##
## Std. Errors:
## (Intercept) Breed_CategoryHerding Group Breed_CategoryHound Group
## Linear 0.2521221 1.109123e+00 0.7341393
## Mixed 70.7807068 7.745849e-04 70.7895356
## Pointed 0.2750003 NaN 1.0957038
## Sharp 70.7807054 7.079248e+01 0.5704836
## Breed_CategoryMixed Breed Breed_CategoryNon-sporting Group
## Linear 6.916461e-01 0.7629374
## Mixed 1.726927e-08 70.7881636
## Pointed 8.370370e-01 0.6578475
## Sharp 1.648954e-03 70.7881615
## Breed_CategorySporting Group Breed_CategoryTerrier Group
## Linear 0.402903219 0.4884307
## Mixed 70.782996110 70.7842389
## Pointed 0.459392449 0.4683738
## Sharp 0.005441496 70.7824715
## Breed_CategoryToy Group Breed_CategoryWorking Group
## Linear 5.746593e-01 0.975942
## Mixed 1.708908e-08 70.788359
## Pointed 5.651339e-01 1.037127
## Sharp 7.078777e+01 777.911696
##
## Residual Deviance: 400.5874
## AIC: 472.5874coefficients3 <- coef(NomModel3)
std_errors3 <- summary_model3$standard.errors
z_values3 <- coefficients3/std_errors3
p_values3 <- 2*(1 - pnorm(abs(z_values3)))
print(p_values3)## (Intercept) Breed_CategoryHerding Group Breed_CategoryHound Group
## Linear 0.0000000 0.0000000 0.0000000
## Mixed 0.7652072 0.0000000 0.7577792
## Pointed 0.0000000 NaN 0.0000000
## Sharp 0.7726918 0.7727269 0.0000000
## Breed_CategoryMixed Breed Breed_CategoryNon-sporting Group
## Linear 0 0.0000000
## Mixed 0 0.7534150
## Pointed 0 0.0000000
## Sharp 0 0.7608653
## Breed_CategorySporting Group Breed_CategoryTerrier Group
## Linear 0.000000 0.0000000
## Mixed 0.745973 0.7503180
## Pointed 0.000000 0.0000000
## Sharp 0.000000 0.7676205
## Breed_CategoryToy Group Breed_CategoryWorking Group
## Linear 0.0000000 0.0000000
## Mixed 0.0000000 0.7727149
## Pointed 0.0000000 0.0000000
## Sharp 0.7553732 0.9652070odds_ratio3 <- exp(summary(NomModel3)$coefficients)## Warning in sqrt(diag(vc)): NaNs produceddata.frame(t(odds_ratio3))## Linear Mixed Pointed
## (Intercept) 5.153714e-10 1.514951e+09 2.356494e-09
## Breed_CategoryHerding Group 3.880018e+09 1.093973e-20 2.650409e-17
## Breed_CategoryHound Group 8.733579e+09 3.301826e-10 2.122673e+08
## Breed_CategoryMixed Breed 3.233827e+09 1.791148e-27 2.828901e+08
## Breed_CategoryNon-sporting Group 1.940334e+09 2.200066e-10 9.901342e+08
## Breed_CategorySporting Group 1.509155e+09 1.100132e-10 1.886046e+08
## Breed_CategoryTerrier Group 2.425388e+09 1.650187e-10 7.956583e+08
## Breed_CategoryToy Group 2.716509e+09 4.639697e-27 6.789790e+08
## Breed_CategoryWorking Group 1.358254e+10 1.320201e-09 1.273095e+09
## Sharp
## (Intercept) 7.574056e+08
## Breed_CategoryHerding Group 1.320094e-09
## Breed_CategoryHound Group 7.996399e-18
## Breed_CategoryMixed Breed 1.552878e-20
## Breed_CategoryNon-sporting Group 4.401092e-10
## Breed_CategorySporting Group 8.540489e-21
## Breed_CategoryTerrier Group 8.251747e-10
## Breed_CategoryToy Group 2.640714e-10
## Breed_CategoryWorking Group 1.832945e-15NomModel4 <- multinom(FB ~ Weight_Category, data = cases)## # weights: 25 (16 variable)
## initial value 259.119504
## iter 10 value 222.046613
## iter 20 value 221.408290
## iter 30 value 221.346929
## final value 221.346500
## convergedsummary_model4 <- summary(NomModel4)
print(summary_model4)## Call:
## multinom(formula = FB ~ Weight_Category, data = cases)
##
## Coefficients:
## (Intercept) Weight_Category> 20 kg Weight_Category10-20 kg
## Linear 0.4044694 -0.2709717 -0.1632883
## Mixed -11.0210631 9.9225312 10.0094901
## Pointed 0.4045527 -0.6158200 -0.3174572
## Sharp -0.6958313 -1.6551891 -1.0086814
## Weight_Category5-10 kg
## Linear 0.3578019
## Mixed -5.7766332
## Pointed 0.1346734
## Sharp 0.1364066
##
## Std. Errors:
## (Intercept) Weight_Category> 20 kg Weight_Category10-20 kg
## Linear 0.9126281 0.9603013 0.9976138
## Mixed 173.7731219 173.7736698 173.7741026
## Pointed 0.9126129 0.9691771 1.0035463
## Sharp 1.2252715 1.4313424 1.4464119
## Weight_Category5-10 kg
## Linear 1.020997
## Mixed 5.771758
## Pointed 1.029109
## Sharp 1.376284
##
## Residual Deviance: 442.693
## AIC: 474.693coefficients4 <- coef(NomModel4)
std_errors4 <- summary_model4$standard.errors
z_values4 <- coefficients4/std_errors4
p_values4 <- 2*(1 - pnorm(abs(z_values4)))
print(p_values4)## (Intercept) Weight_Category> 20 kg Weight_Category10-20 kg
## Linear 0.6576269 0.7778104 0.8699839
## Mixed 0.9494304 0.9544653 0.9540668
## Pointed 0.6575556 0.5251642 0.7517480
## Sharp 0.5701031 0.2475219 0.4855725
## Weight_Category5-10 kg
## Linear 0.7260059
## Mixed 0.3169019
## Pointed 0.8958829
## Sharp 0.9210491odds_ratio4 <- exp(summary(NomModel4)$coefficients)
data.frame(t(odds_ratio4))## Linear Mixed Pointed Sharp
## (Intercept) 1.4985073 1.635359e-05 1.4986320 0.4986597
## Weight_Category> 20 kg 0.7626380 2.038452e+04 0.5401977 0.1910559
## Weight_Category10-20 kg 0.8493463 2.223649e+04 0.7279978 0.3646995
## Weight_Category5-10 kg 1.4301822 3.099132e-03 1.1441630 1.1461479NomModel4 <- multinom(FB ~ Breed_Category + Age_Category + Weight_Category + Sex + Neuter, data = cases)## # weights: 90 (68 variable)
## initial value 259.119504
## iter 10 value 191.793358
## iter 20 value 181.736397
## iter 30 value 179.280036
## iter 40 value 178.765448
## iter 50 value 178.687744
## iter 60 value 178.684077
## final value 178.684036
## convergedsummary_model4 <- summary(NomModel4)
print(summary_model4)## Call:
## multinom(formula = FB ~ Breed_Category + Age_Category + Weight_Category +
## Sex + Neuter, data = cases)
##
## Coefficients:
## (Intercept) Breed_CategoryHerding Group Breed_CategoryHound Group
## Linear -12.331631 12.73565 13.65795
## Mixed -8.729909 -29.62492 -13.35493
## Pointed -12.532079 -17.92785 10.61710
## Sharp -10.491893 -14.04057 -38.10596
## Breed_CategoryMixed Breed Breed_CategoryNon-sporting Group
## Linear 12.45687 11.64172
## Mixed -30.96014 -13.72473
## Pointed 10.76086 12.21848
## Sharp -36.67760 -15.46548
## Breed_CategorySporting Group Breed_CategoryTerrier Group
## Linear 11.72947 12.03081
## Mixed -15.40045 -14.40971
## Pointed 10.58606 11.75612
## Sharp -34.60294 -15.19740
## Breed_CategoryToy Group Breed_CategoryWorking Group Age_Category>= 8 Y
## Linear 12.10489 13.87592 0.7968135
## Mixed -28.12534 -12.97112 0.3848321
## Pointed 11.56708 12.47456 1.9725120
## Sharp -17.22661 -30.07859 0.2918869
## Age_Category2-5 Y Age_Category5-8 Y Weight_Category> 20 kg
## Linear -0.2273787 0.5665018 0.6500772
## Mixed -1.4011911 -22.5495879 23.1791236
## Pointed 0.8467703 1.3369056 -0.7073701
## Sharp 0.8098021 0.9317736 -2.1342947
## Weight_Category10-20 kg Weight_Category5-10 kg SexM NeuterY
## Linear 0.7473795 1.1729831 -0.3658228 -0.3336033
## Mixed 22.9750998 -10.2593899 0.0441858 -0.7945208
## Pointed -0.5808108 -0.1354941 0.1476527 0.8632696
## Sharp -2.1761921 0.3271637 -1.8193824 26.6371899
##
## Std. Errors:
## (Intercept) Breed_CategoryHerding Group Breed_CategoryHound Group
## Linear 0.9797944 1.193857e+00 7.687640e-01
## Mixed 0.7989772 5.520773e-07 1.508517e+00
## Pointed 1.0745543 7.996800e-13 1.133341e+00
## Sharp 1.0044280 1.676463e+00 3.611747e-09
## Breed_CategoryMixed Breed Breed_CategoryNon-sporting Group
## Linear 7.729064e-01 0.8153067
## Mixed 1.438647e-07 1.4188557
## Pointed 9.148957e-01 0.7198376
## Sharp 8.297823e-09 1.4417883
## Breed_CategorySporting Group Breed_CategoryTerrier Group
## Linear 5.397865e-01 0.5371527
## Mixed 1.083729e+00 1.1680906
## Pointed 6.268750e-01 0.5424714
## Sharp 5.905645e-08 1.2724736
## Breed_CategoryToy Group Breed_CategoryWorking Group Age_Category>= 8 Y
## Linear 6.893544e-01 1.019262e+00 0.6961466
## Mixed 1.082871e-06 1.478039e+00 1.0620854
## Pointed 6.772895e-01 1.130259e+00 0.7568055
## Sharp 1.790926e+00 4.727599e-07 1.4150059
## Age_Category2-5 Y Age_Category5-8 Y Weight_Category> 20 kg
## Linear 0.5881065 7.893121e-01 1.241854
## Mixed 1.2604365 6.357843e-10 0.631826
## Pointed 0.6842283 8.707181e-01 1.285559
## Sharp 1.3309494 1.528639e+00 1.903266
## Weight_Category10-20 kg Weight_Category5-10 kg SexM NeuterY
## Linear 1.2087801 1.181085e+00 0.4996435 0.5961217
## Mixed 0.5856722 3.307195e-14 0.9279511 0.9269251
## Pointed 1.2540538 1.202356e+00 0.5302742 0.7684698
## Sharp 1.9927784 1.729339e+00 1.0711851 1.0044280
##
## Residual Deviance: 357.3681
## AIC: 493.3681coefficients4 <- coef(NomModel4)
std_errors4 <- summary_model4$standard.errors
z_values4 <- coefficients4/std_errors4
p_values4 <- 2*(1 - pnorm(abs(z_values4)))
print(p_values4)## (Intercept) Breed_CategoryHerding Group Breed_CategoryHound Group
## Linear 0 0 0
## Mixed 0 0 0
## Pointed 0 0 0
## Sharp 0 0 0
## Breed_CategoryMixed Breed Breed_CategoryNon-sporting Group
## Linear 0 0
## Mixed 0 0
## Pointed 0 0
## Sharp 0 0
## Breed_CategorySporting Group Breed_CategoryTerrier Group
## Linear 0 0
## Mixed 0 0
## Pointed 0 0
## Sharp 0 0
## Breed_CategoryToy Group Breed_CategoryWorking Group Age_Category>= 8 Y
## Linear 0 0 0.252372473
## Mixed 0 0 0.717100687
## Pointed 0 0 0.009150864
## Sharp 0 0 0.836572477
## Age_Category2-5 Y Age_Category5-8 Y Weight_Category> 20 kg
## Linear 0.6990313 0.4729325 0.6006451
## Mixed 0.2662795 0.0000000 0.0000000
## Pointed 0.2158810 0.1246842 0.5821525
## Sharp 0.5428961 0.5421635 0.2621239
## Weight_Category10-20 kg Weight_Category5-10 kg SexM NeuterY
## Linear 0.5363826 0.3206415 0.46406630 0.5757367
## Mixed 0.0000000 0.0000000 0.96202187 0.3913578
## Pointed 0.6432593 0.9102759 0.78067003 0.2612839
## Sharp 0.2748159 0.8499484 0.08941793 0.0000000joint_tests(NomModel4)## Error in solve.default(zcov, z) :
## system is computationally singular: reciprocal condition number = 5.9088e-19
## Error in solve.default(zcov, z) :
## system is computationally singular: reciprocal condition number = 1.27874e-23
## Error in solve.default(zcov, z) :
## system is computationally singular: reciprocal condition number = 1.87768e-20## model term df1 df2 F.ratio p.value note
## Breed_Category 8 68 0.000 1.0000
## Age_Category 3 68 0.000 1.0000
## Weight_Category 3 68 0.000 1.0000
## Sex 1 68 0.000 1.0000
## Neuter 1 68 0.000 1.0000
## FB 4 68 35.461 <.0001
## Breed_Category:FB 32 NA NA NA d
## Age_Category:FB 12 68 2.866 0.0030
## Weight_Category:FB 12 68 2.266 0.0175
## Sex:FB 4 68 1.317 0.2724
## Neuter:FB 4 68 7.164 0.0001
## (confounded) 2827 NA NA NA
##
## d: df1 reduced due to linear dependenceodds_ratio4 <- exp(summary(NomModel4)$coefficients)
data.frame(t(odds_ratio4))## Linear Mixed Pointed
## (Intercept) 4.410022e-06 1.616772e-04 3.609004e-06
## Breed_CategoryHerding Group 3.396433e+05 1.361640e-13 1.636936e-08
## Breed_CategoryHound Group 8.542282e+05 1.585000e-06 4.082713e+04
## Breed_CategoryMixed Breed 2.570091e+05 3.582455e-14 4.713902e+04
## Breed_CategoryNon-sporting Group 1.137458e+05 1.095028e-06 2.024965e+05
## Breed_CategorySporting Group 1.241778e+05 2.049610e-07 3.957920e+04
## Breed_CategoryTerrier Group 1.678475e+05 5.520017e-07 1.275311e+05
## Breed_CategoryToy Group 1.807538e+05 6.099852e-13 1.055643e+05
## Breed_CategoryWorking Group 1.062276e+06 2.326568e-06 2.615979e+05
## Age_Category>= 8 Y 2.218461e+00 1.469368e+00 7.188712e+00
## Age_Category2-5 Y 7.966191e-01 2.463034e-01 2.332103e+00
## Age_Category5-8 Y 1.762092e+00 1.610046e-10 3.807244e+00
## Weight_Category> 20 kg 1.915689e+00 1.165643e+10 4.929389e-01
## Weight_Category10-20 kg 2.111460e+00 9.505152e+09 5.594446e-01
## Weight_Category5-10 kg 3.231619e+00 3.502705e-05 8.732843e-01
## SexM 6.936257e-01 1.045177e+00 1.159110e+00
## NeuterY 7.163379e-01 4.517977e-01 2.370900e+00
## Sharp
## (Intercept) 2.776060e-05
## Breed_CategoryHerding Group 7.984648e-07
## Breed_CategoryHound Group 2.823523e-17
## Breed_CategoryMixed Breed 1.177932e-16
## Breed_CategoryNon-sporting Group 1.920565e-07
## Breed_CategorySporting Group 9.378476e-16
## Breed_CategoryTerrier Group 2.511034e-07
## Breed_CategoryToy Group 3.300497e-08
## Breed_CategoryWorking Group 8.650390e-14
## Age_Category>= 8 Y 1.338952e+00
## Age_Category2-5 Y 2.247463e+00
## Age_Category5-8 Y 2.539008e+00
## Weight_Category> 20 kg 1.183280e-01
## Weight_Category10-20 kg 1.134728e-01
## Weight_Category5-10 kg 1.387029e+00
## SexM 1.621258e-01
## NeuterY 3.701558e+11