Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit


Answers for MA6521 2021 Exam

Note there is often more than one possible answer to each question. The answers are illustrative examples. Note possible to use either SPSS or R as applicable.

1a)

The length of stay in hospital - complexity of operation.

Firstly test for normality:

H0: Length of stay in hospital (split by complex) follows a normal distribution.

H1: Length of stay in hospital (split by complex) does not follow a normal distribution.

Evidence to reject null hypothesis in some cases, but not others. Better to be cautious here and assume data does not follow normal distribution.

Then we test whether there is a difference of length of stay for simple and complex operations:       H0: There is no difference between length of stay in hospital regardless of complexity of operation. H1: There is a difference between length of stay in hospital depending on complexity of operation.

There is strong evidence to reject the null hypothesis (p-value <0.001).

 


On average the length of stay in hospital is significantly larger ifthe operation was more complex.

There is a large difference here, so the sample size is large enough to show there is a significant difference.

The length of stay in hospital - BMI

Check whether follows a normal distribution.

H0: Length of stay in hospital (split by BMI) follows a normal distribution.

H1: Length of stay in hospital (split by BMI) does not follow a normal distribution.

Do not reject the null hypothesis, a parametric test would be sufficient in this case.  However, could also justify using a non-parametric test as data set is too small and not enough power to pick up data is non-normal. Marks will be        given based on sensible discussion that is justified.

H0: There is no difference in length of stay in hospital split by BMI

H1: There is a difference in length of stay in hospital split by BMI

Or

 

Do not reject the null hypothesis in either case.

 

 


 

There is a small difference in mean length of stay in hospital with BMI category 1 having lowest mean (4.33), then      category 2 (mean 5.14), then category 3 (mean 5.88) . However, there is not a significant difference.  A larger sample size would be needed to show whether this difference exists or is due to random change.

A power analysis is needed to estimate sample size. Students would need to assume data followed a normal distribution to do this.

R code:

> library(pwr)

> cohen.ES(test="anov",size="small")

Conventional effect size from Cohen (1982)

test = anov


size = small

effect.size = 0.1

> pwr.anova.test(k =3 , f = 0.1 , sig.level = 0.05, power = 0.8)

Balanced one-way analysis of variance power calculation

k = 3

n = 322.157

f = 0.1

sig.level = 0.05

power = 0.8

NOTE: n is number in each group

> cohen.ES(test="anov",size="medium")

Conventional effect size from Cohen (1982)

test = anov

size = medium

effect.size = 0.25

> pwr.anova.test(k =3 , f = 0.25 , sig.level = 0.05, power = 0.8)

Balanced one-way analysis of variance power calculation

k = 3

n = 52.3966

f = 0.25

sig.level = 0.05

power = 0.8

NOTE: n is number in each group

> pp <- c(sum(data1$BMI ==1),sum(data1$BMI ==2),sum(data1$BMI ==3)) > p <- pp/sum(pp)

> gmu <- c(mean(data1$LOS[data1$BMI==1]),mean(data1$LOS[data1$BMI==2]), mean(data1$LOS[data1$BMI==3]))

> mu <- mean(data1$LOS)

> ff <- sqrt(sum(p*(gmu-mu)^2))/sd(data1$LOS)

> pwr.anova.test(k =3 , f = ff , sig.level = 0.05, power = 0.8)

Balanced one-way analysis of variance power calculation

k = 3

n = 59.11045

f = 0.2351088

sig.level = 0.05

power = 0.8

NOTE: n is number in each group

 

A small effect size needs a sample size of 323 per group, a medium 52.4 per group. Or if we had the exact mean differences from the sample a sample size of 60 in each group is need, so 180 in total.


b) This report is about data on length of stay in hospital.

There is strong evidence in the data that a more complex operation leads to a longer stay in hospital. The average stay in hospital for a simple operation is 3.1 days whereas it is 7.0 days for a complex operation. This is further       illustrated in the figure below, which shows complex operations have longer stays than simple operations.

 

Length of Stay in Hospital by Complexity of Operation

40%

35%

30%

25%

20%

15%

10%

5%

0%

1         2         3         4         5         6         7         8         9        10

Length of Stay, in Days

 

 

There is no evidence in the data that a patient’s BMI effects their length of stay in hospital. Whilst the mean for BMI less than 25.0 is 4.33, for BMI 25.0 to 30.0 the mean is 5.14 and for BMI greater or equal to 3.00 the mean is 5.88,     there is no evidence in the data that this is not just random chance. A larger sample size of 60 patients in each BMI   category, 180 patients in total, would be needed to detect a difference of this size.

In conclusion the length of stay in hospital does depend on how complex an operation is, but not the BMI.


2) Fit a Poisson GLM as number of elephants is a count.

Mit is a factor, temp is a continuous covariate and loc is a random effect.

First consider all variables except random effect:

R code to open data set:

> library(readxl)

> datael <- read_excel("exam_MA6512_2021.xlsx",sheet=2)

R code to fit the full model (without random effect).

R code to fit the other combinations of models

 

> modelm<-glm(nel~ > modelt<-glm(nel~ > modelc<-glm(nel~

factor(mit),data=datael,family = poisson()) temp,data=datael,family = poisson())        1,data=datael,family = poisson())

 

Model selection GLMs AIC:

> AICvalues <- c(AIC(modelmt),AIC(modelm),AIC(modelt),AIC(modelc)) > modelname <- c("mit+temp", "mit", "temp","const")

> AICvaluesr <- round(AICvalues,digits=2)

> deltaAIC <-AICvalues-min(AICvalues)

> AICweight <- round(exp(-deltaAIC/2)/sum(exp(-deltaAIC/2)),digits=2) > deltaAIC <-round(deltaAIC,digits=2)

> AICtable = data.frame(modelname,AICvaluesr,deltaAIC,AICweight)    > colnames(AICtable ) <- c("Model", "AIC","Delta AIC","AIC weight") > print(AICtable)

 


Model 1 mit+temp

2      mit

3     temp

4    const


AIC Delta AIC AIC weight

3477.88

3477.37

5178.07

5181.08


The best model, the one with the smallest AIC, includes mit but not time. The model weight is 0.56, mit+temp is close with a model weight of 0.44.

R code for model selection analysis of deviance.

> # Model selection Deviance

> # Level 1 - compare to constant model

> Modellist <- list( )

> Modellist[[1]] <- modelm

> Modellist[[2]] <- modelt

> modelcompare <- modelc

> modelname <- c("mit", "temp")

> dev <- rep(0,(length(Modellist)))

> for (i in 1:(length(Modellist))){dev[i]=Modellist[[i]]$deviance} > devcompare <- modelcompare$deviance

> DF <- rep(0,(length(Modellist)))

> for (i in 1:(length(Modellist))){DF[i]=Modellist[[i]]$df.residual} > DFcompare <- modelcompare$df.residual

> diffdev=devcompare-dev

> diffDF=DFcompare-DF

> pval=pchisq(diffdev, df=diffDF, lower.tail=FALSE)

> devtable = data.frame(modelname,diffdev,diffDF,pval)

> colnames(devtable) <- c("Model", "Diff Deviance","Diff DF","p-value") > print(devtable)

 


Model Diff Deviance Diff DF

1709.713092

5.010123


p-value

0.00000000

0.02519951


At level 1 (comparing to constant model) both make a significant difference in deviance, mit has smallest p-value so compare at next stage.

At level 2 (comparing to mit) temp is not significant. Therefore the best model includes only mit.             Both AIC and analysis of deviance result in the same conclusion that best model has only the factor mit

R code for fitting mit + location random effect

> library(lme4)

> modelmr <- glmer(nel~ factor(mit)+(1|loc),data=datael,family = poisson())

LRT testing whether need random effect

> LRT <- 2*(logLik(modelmr)[1]-logLik(modelm)[1])

> pvalue <- 1-pchisq(LRT, df=1)


> LRTresults <- data.frame(LRT,pvalue)

> print(LRTresults)

LRT pvalue

1 117.5948      0

 

AIC model selection

> library(lmerTest)

> AIC(modelmr,modelm)

df      AIC

modelmr  5 3361.771

modelm   4 3477.366

Both AIC and LRT support a random effect

Output for best model

> summary(modelmr)

Generalized linear mixed model fit by maximum likelihood (Laplace

Approximation) [glmerMod]

Family: poisson  ( log )

Formula: nel ~ factor(mit) + (1 | loc)

Data: datael

 

 


AIC 3361.8


BIC   logLik deviance df.resid

3386.5  -1675.9   3351.8     1025


 

Scaled residuals:

Min      1Q  Median      3Q     Max

-1.9242 -0.5998 -0.2972  0.5834  4.9739

Random effects:

Groups Name        Variance Std.Dev.

loc    (Intercept) 0.05312  0.2305

Number of obs: 1030, groups:  obs, 15

Fixed effects:

Estimate Std. Error z value Pr(>|z|)

 


(Intercept)  factor(mit)1 factor(mit)2 factor(mit)3 ---


1.77049 -2.86108 -0.71153 -1.27037


0.06496 0.10810 0.04411 0.05728


27.25 -26.47 -16.13 -22.18


<2e-16 ***

<2e-16 ***

<2e-16 ***

<2e-16 ***


Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:

(Intr) fct()1 fct()2

 


factor(mt)1 factor(mt)2 factor(mt)3


-0.096 -0.236 -0.181


0.144

0.111  0.271


 

The best model includes the mitigation and the random location effect .

The parameter estimates and standard errors for this model are as follows:

Parameter                   Estimate        SE            

Intercept                     1.77               0.07

Mit 1                            -2.86              0.11

Mit 2                            -0.71              0.04

Mit 3                            -1.27              0.06

Location RE                0.23

(variance)

The number of elephants seen depends on the mitigation used. No mitigation (0) is when the most elephants are       seen followed by mitigation 2, then mitigation 3 with mitigation 1 when the least elephants are seen. The number of elephants also varies with the random location effect.