Flight Landing Distance Analysis
For binary response using Logistic Regression
Vidhi Rathod
13 February 2020
Landing Assessment
1.0 Case Overview:
From now on, we work on the cleaned FAA data set prepared by carrying out the prior project using Linear regression. For other readers I am including the data load part again:
Here,We use same FAA data to predict two types of landings: Risky and long
Factors Affecting Actual Landing Distance
Variable dictionary:
- Aircraft: The make of an aircraft (Boeing or Airbus).
- Duration (in minutes): Flight duration between taking off and landing. The duration of a normal flight should always be greater than 40min.
- No_pasg: The number of passengers in a flight.
- Speed_ground (in miles per hour): The ground speed of an aircraft when passing over the threshold of the runway. If its value is less than 30MPH or greater than 140MPH, then the landing would be considered as abnormal.
- Speed_air (in miles per hour): The air speed of an aircraft when passing over the threshold of the runway. If its value is less than 30MPH or greater than 140MPH, then the landing would be considered as abnormal.
- Height (in meters): The height of an aircraft when it is passing over the threshold of the runway. The landing aircraft is required to be at least 6 meters high at the threshold of the runway.
- Pitch (in degrees): Pitch angle of an aircraft when it is passing over the threshold of the runway.
- Distance (in feet): The landing distance of an aircraft. More specifically, it refers to the distance between the threshold of the runway and the point where the aircraft can be fully stopped. The length of the airport runway is typically less than 6000 feet.
2.1 Initial Set Up
2.2 Package Loading
Required Packages
- Tidyverse (dplyr, ggplot2..) - Data Read, Manipulation and visualisation
- Plotly - Interactive Visualization
- KableExtra - Styling for table (Styling Data Tables within Markdown)
- gridExtra- Graphical arrangement
- forecast- Time series and forecasting
- ggplot2- Graphical representation
- stringr- string manipulations
- corrplot-making correlogram
- knitrr- Dynamic report generation
2.3 Data Loading
Data Load
# Read the sheets, one by one
FAA1 <- read_excel("FAA1.xls")
FAA2 <- read_excel("FAA2.xls")
dim(FAA1)## [1] 800 8dim(FAA2)## [1] 150 7Structure
str(FAA1)## Classes 'tbl_df', 'tbl' and 'data.frame': 800 obs. of 8 variables:
## $ aircraft : chr "boeing" "boeing" "boeing" "boeing" ...
## $ duration : num 98.5 125.7 112 196.8 90.1 ...
## $ no_pasg : num 53 69 61 56 70 55 54 57 61 56 ...
## $ speed_ground: num 107.9 101.7 71.1 85.8 59.9 ...
## $ speed_air : num 109 103 NA NA NA ...
## $ height : num 27.4 27.8 18.6 30.7 32.4 ...
## $ pitch : num 4.04 4.12 4.43 3.88 4.03 ...
## $ distance : num 3370 2988 1145 1664 1050 ...str(FAA2)## Classes 'tbl_df', 'tbl' and 'data.frame': 150 obs. of 7 variables:
## $ aircraft : chr "boeing" "boeing" "boeing" "boeing" ...
## $ no_pasg : num 53 69 61 56 70 55 54 57 61 56 ...
## $ speed_ground: num 107.9 101.7 71.1 85.8 59.9 ...
## $ speed_air : num 109 103 NA NA NA ...
## $ height : num 27.4 27.8 18.6 30.7 32.4 ...
## $ pitch : num 4.04 4.12 4.43 3.88 4.03 ...
## $ distance : num 3370 2988 1145 1664 1050 ...FAA1 : 800 Rows and 8 columns , sample size is 800 rows
FAA2 : 150 Rows and 7 columns ,sample size is 150 rows
FAA2 does not have duration column
FAA1 AND FAA2 have similar structure, datatypes.
4.4 Initial Data processing
In order to combine the 2 data sets we need to make the sructure of both same, by structure we refer to columns.Hence we need to add a column duration to FAA2 and then use rbind function to find 150 +800 =950 rows sample size=950
Data Merge
Merge the 2 data sets FAA1 & 2 to create a composite dataset
FAA2$duration<-NA
FAA_merge<-rbind(FAA1,FAA2)Duplicate Check: Because the FAA2 has NA in duration column , we need to check for duplicates rows excluding that column. We need to remove the duplicate rows if we find any to reduce redundancy
row_dup<-duplicated(FAA_merge[,-2])
class(row_dup)## [1] "logical"FAA_U<-FAA_merge[!row_dup,]
dim(FAA_U)## [1] 850 8Summary Statistics
we find 850 unique rows and 8 columns. Below is the summary statistics of each column in dataset
dim(FAA_U)## [1] 850 8str(FAA_U)## Classes 'tbl_df', 'tbl' and 'data.frame': 850 obs. of 8 variables:
## $ aircraft : chr "boeing" "boeing" "boeing" "boeing" ...
## $ duration : num 98.5 125.7 112 196.8 90.1 ...
## $ no_pasg : num 53 69 61 56 70 55 54 57 61 56 ...
## $ speed_ground: num 107.9 101.7 71.1 85.8 59.9 ...
## $ speed_air : num 109 103 NA NA NA ...
## $ height : num 27.4 27.8 18.6 30.7 32.4 ...
## $ pitch : num 4.04 4.12 4.43 3.88 4.03 ...
## $ distance : num 3370 2988 1145 1664 1050 ...summary(FAA_U)## aircraft duration no_pasg speed_ground
## Length:850 Min. : 14.76 Min. :29.0 Min. : 27.74
## Class :character 1st Qu.:119.49 1st Qu.:55.0 1st Qu.: 65.90
## Mode :character Median :153.95 Median :60.0 Median : 79.64
## Mean :154.01 Mean :60.1 Mean : 79.45
## 3rd Qu.:188.91 3rd Qu.:65.0 3rd Qu.: 92.06
## Max. :305.62 Max. :87.0 Max. :141.22
## NA's :50
## speed_air height pitch distance
## Min. : 90.00 Min. :-3.546 Min. :2.284 Min. : 34.08
## 1st Qu.: 96.25 1st Qu.:23.314 1st Qu.:3.642 1st Qu.: 883.79
## Median :101.15 Median :30.093 Median :4.008 Median :1258.09
## Mean :103.80 Mean :30.144 Mean :4.009 Mean :1526.02
## 3rd Qu.:109.40 3rd Qu.:36.993 3rd Qu.:4.377 3rd Qu.:1936.95
## Max. :141.72 Max. :59.946 Max. :5.927 Max. :6533.05
## NA's :642INSIGHTS:
- Duration column in FAA2 file is missing
- Speed air column has 642 NA Values
- duration,distance and height columns have considerable difference between Minimum and Maximum value
- height column has a negative value, which cannot be legit because height is bound to be positive.Hence it can be a data reading issue.
3.0 Data Cleaning
3.1 Initial scrub
As first steps in Data cleaning we need to filter the rows that do not qualify for data analysis.The are 2 major categories of such line items: * Abnormal values: We need to remove abnormal values(values that do not qualify for analysis, as suggested by SMEs).Below code removes them * NA Rows: We need to remove NA rows(where all cell values are NA).
FAA_U<-FAA_U[FAA_U$duration>=40,]
FAA_U<-FAA_U[FAA_U$speed_ground>=30 && FAA_U$speed_ground<=140,]
FAA_U<-FAA_U[FAA_U$height>=6,]
FAA_U<-FAA_U[FAA_U$distance<=6000,]
nrow(FAA_U)## [1] 833sum(is.na(FAA_U$duration))## [1] 50sum(is.na(FAA_U$duration))## [1] 50sum(is.na(FAA_U$speed_ground))## [1] 50sum(is.na(FAA_U$height))## [1] 50sum(is.na(FAA_U$distance))## [1] 50#Remove all rowa that have NA in all the columns
FAA_U_RmvNA<-FAA_U %>% filter(!(is.na(aircraft)&is.na(duration)&is.na(no_pasg)&is.na(speed_ground)&is.na(speed_air)&is.na(height)&is.na(pitch)&is.na(distance)))
#Looking for NA Values in columns
colSums(is.na(FAA_U_RmvNA))## aircraft duration no_pasg speed_ground speed_air
## 0 0 0 0 588
## height pitch distance
## 0 0 0- Conclusion: we removed close to 850 -783= 67 rows, ** Abnormal values: 17 ** NA VALUE Rows: 50
3.2 Binary Response
Modeling a binary response using logistic regression Step 01
Create two binary variables below and attach them to your data set.
long.landing = 1 if distance > 2500; =0 otherwise risky.landing = 1 if distance > 3000; =0 otherwise.
Discard the continuous data you have for “distance”, and assume we are given the binary data of “long.landing” and “risky.landing” only.
Use above formula to create 2 new columns
FAA_U_RmvNA$long.landing<-ifelse(FAA_U_RmvNA$distance>=2500,1,0)
FAA_U_RmvNA$risky.landing<-ifelse(FAA_U_RmvNA$distance>=3000,1,0)We will not be using the distance variable going forward
4.0 Data Analysis
4.1 Long Landing
Univariate analysis
STEP2
Identifying important factors using the binary data of “long.landing”.
Use a pie chart or a histogram to show the distribution of “long.landing”.
par(mfrow=c(2,1))
hist(FAA_U_RmvNA$long.landing,main = paste("Histogram of long.landing"),xlab="long.landing")
#Univariate analysis:risky.landing
hist(FAA_U_RmvNA$long.landing,main = paste("Histogram of long.landing"),xlab="risky.landing")
Step 3 * Perform single-factor regression analysis for each of the potential risk factors, in a similar way to what you did in Steps 13-15 of Part 1. But here the response “long.landing” is binary. You may consider using logistic regression. Provide a table that ranks the factors from the most important to the least. This table contains 5 columns: the names of variables, the size of the regression coefficient, the odds ratio, the direction of the regression coefficient (positive or negative), and the p-value.
#names(FAA_U_RmvNA)
FAAlog<-FAA_U_RmvNA[,-8]
### A model with all variables available
lmod.aircraft<-glm(long.landing~aircraft,family=binomial,FAAlog)
lmod.duration<-glm(long.landing~duration,family=binomial,FAAlog)
lmod.no_pasg<-glm(long.landing~no_pasg,family=binomial,FAAlog)
lmod.speed_ground<-glm(long.landing~speed_ground,family=binomial,FAAlog)
lmod.speed_air<-glm(long.landing~speed_air,family=binomial,FAAlog)
lmod.height<-glm(long.landing~height,family=binomial,FAAlog)
lmod.pitch<-glm(long.landing~pitch,family=binomial,FAAlog)
v1<-c('aircraft',round(coef(lmod.aircraft)[1],2),coef(lmod.aircraft)[2],exp(coef(lmod.aircraft))[1],exp(coef(lmod.aircraft))[2],
sign<-ifelse(coef(lmod.aircraft)[2]==-1,"Negative","Positive"),coef(summary(lmod.aircraft))[,4])
v2<-c('duration',coef(lmod.duration)[1],coef(lmod.duration)[2],exp(coef(lmod.duration))[1],exp(coef(lmod.duration))[2],
sign<-ifelse(coef(lmod.duration)[2]==-1,"Negative","Positive"),coef(summary(lmod.duration))[,4])
v3<-c('no_pasg',coef(lmod.no_pasg)[1],coef(lmod.no_pasg)[2],exp(coef(lmod.no_pasg))[1],exp(coef(lmod.no_pasg))[2],
sign<-ifelse(sign(coef(lmod.no_pasg)[2])==-1,"Negative","Positive"),coef(summary(lmod.no_pasg))[,4])
v4<-c('speed_ground',coef(lmod.speed_ground)[1],coef(lmod.speed_ground)[2],exp(coef(lmod.speed_ground))[1],exp(coef(lmod.speed_ground))[2],
sign<-ifelse(coef(lmod.speed_ground)[2]==-1,"Negative","Positive"),coef(summary(lmod.speed_ground))[,4])
v5<-c('speed_air',coef(lmod.speed_air)[1],coef(lmod.speed_air)[2],exp(coef(lmod.speed_air))[1],exp(coef(lmod.speed_air))[2],
sign<-ifelse(coef(lmod.speed_air)[2]==-1,"Negative","Positive"),coef(summary(lmod.speed_air))[,4])
v6<-c('height',coef(lmod.height)[1],coef(lmod.height)[2],exp(coef(lmod.height))[1],exp(coef(lmod.height))[2],
sign<-ifelse(coef(lmod.height)[2]==-1,"Negative","Positive"),coef(summary(lmod.height))[,4])
v7<-c('pitch',coef(lmod.pitch)[1],coef(lmod.pitch)[2],exp(coef(lmod.pitch))[1],exp(coef(lmod.pitch))[2],
sign<-ifelse(coef(lmod.pitch)[2]==-1,"Negative","Positive"),coef(summary(lmod.pitch))[,4])
d<-data.frame(rbind(v1,v2,v3,v4,v5,v6,v7))
colnames(d)<-c('Name','B0','B1','OddsRto_B0','OddsRto_B1','Sign','Pvalue_B0','Pvalue_B1')
d$B0<-round(as.numeric(as.character(d$B0)),2)
d$B1<-round(as.numeric(as.character(d$B1)),2)
d$OddsRto_B0<-round(as.numeric(as.character(d$OddsRto_B0)),2)
d$OddsRto_B1<-round(as.numeric(as.character(d$OddsRto_B1)),2)
d$Pvalue_B0<-as.numeric(as.character(d$Pvalue_B0))
d$Pvalue_B1<-as.numeric(as.character(d$Pvalue_B1))
d<-arrange(d,Pvalue_B1,Pvalue_B0)
d$signif<-ifelse(d$Pvalue_B1<0.05,'S','I')
d[,c(1,2,3,6,7,8,9)]## Name B0 B1 Sign Pvalue_B0 Pvalue_B1 signif
## 1 speed_ground -48.30 0.48 Positive 7.492101e-14 1.296677e-13 S
## 2 speed_air -52.89 0.52 Positive 1.187237e-10 1.514441e-10 S
## 3 aircraft -2.39 0.82 Positive 1.592773e-39 2.743033e-04 S
## 4 pitch -3.25 0.33 Positive 1.277356e-04 1.107887e-01 I
## 5 height -2.14 0.01 Positive 1.846205e-09 5.128363e-01 I
## 6 duration -1.75 0.00 Positive 9.116182e-07 6.213848e-01 I
## 7 no_pasg -1.53 -0.01 Negative 7.471448e-02 6.438361e-01 IObservation:
- 3 predictors seem to be significant based on pvalue analysis, namely: speed_ground,speed_air,aircraft
- All coefficients are positive: which means they positively effect the landing distance
- Aircraft has the biggest coefficient follwed by speed air and speed ground.
- We also observe that pvalue for height is not significan.
Step 4. For those significant factors identified in Step 3, visualize its association with “long.landing”.
FAAlog$aircraft1<-ifelse(FAAlog$aircraft=='boeing',as.numeric(1),as.numeric(0))
par(mfrow=c(1,2))
plot(long.landing~speed_ground,FAAlog)
plot(jitter(long.landing,0.2)~jitter(speed_ground),FAAlog,xlab ="speed ground",ylab ="long.landing",pch=".")
par(mfrow=c(1,2))
plot(long.landing~speed_air,FAAlog)
plot(jitter(long.landing,0.2)~jitter(speed_air),FAAlog,xlab ="speed_air",ylab ="long.landing",pch=".")
par(mfrow=c(1,2))
plot(long.landing~aircraft1,FAAlog)
plot(jitter(long.landing,0.2)~jitter(aircraft1),FAAlog,xlab ="aircraft",ylab ="long.landing",pch=".")
Fit a model
Step 5. Based on the analysis results in Steps 3-4 and the collinearity result seen in earlier, initiate a “full” model. Fit your model to the data and present your result.
Understanding: We observe that speed air and speed ground both are associated with Landing distance, these predictors are have have high magnitude of multi collinearity amongst them. We understand that picking one out of them can help us remove multicollinearity from our model.
Here we choose speed ground rather that speed air because speed air have almost 75% missing values.
#full model
lmod.full_faa<-glm(long.landing~aircraft1+duration+no_pasg+speed_ground+height+pitch,family=binomial,FAAlog)
summary(lmod.full_faa)##
## Call:
## glm(formula = long.landing ~ aircraft1 + duration + no_pasg +
## speed_ground + height + pitch, family = binomial, data = FAAlog)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.12757 -0.00078 0.00000 0.00000 2.19551
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.149e+02 2.411e+01 -4.765 1.89e-06 ***
## aircraft1 4.985e+00 1.181e+00 4.221 2.44e-05 ***
## duration 5.361e-03 7.704e-03 0.696 0.4865
## no_pasg 7.420e-03 5.559e-02 0.133 0.8938
## speed_ground 9.795e-01 2.006e-01 4.883 1.05e-06 ***
## height 2.346e-01 7.188e-02 3.264 0.0011 **
## pitch 1.283e+00 8.423e-01 1.523 0.1278
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 598.240 on 782 degrees of freedom
## Residual deviance: 51.084 on 776 degrees of freedom
## AIC: 65.084
##
## Number of Fisher Scoring iterations: 12Upon regressing full model we realize that aircraft, speed_ground & height are significant as per p values.
Doing a further drop1 analysis can help us be confident about this model:
fit a full model after removing collinerity
drop1 analysis
#######################################################################
#fit a full model after removing collinerity
drop1(lmod.full_faa,test ="Chi")## Single term deletions
##
## Model:
## long.landing ~ aircraft1 + duration + no_pasg + speed_ground +
## height + pitch
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.08 65.08
## aircraft1 1 85.82 97.82 34.73 3.785e-09 ***
## duration 1 51.57 63.57 0.49 0.4836
## no_pasg 1 51.10 63.10 0.02 0.8937
## speed_ground 1 583.37 595.37 532.29 < 2.2e-16 ***
## height 1 72.65 84.65 21.57 3.414e-06 ***
## pitch 1 53.67 65.67 2.59 0.1076
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(lmod.full_faa, ~ . -no_pasg ), test = "Chi")## Single term deletions
##
## Model:
## long.landing ~ aircraft1 + duration + speed_ground + height +
## pitch
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.10 63.10
## aircraft1 1 85.86 95.86 34.76 3.724e-09 ***
## duration 1 51.58 61.58 0.48 0.4894
## speed_ground 1 583.54 593.54 532.44 < 2.2e-16 ***
## height 1 73.47 83.47 22.36 2.257e-06 ***
## pitch 1 53.68 63.68 2.58 0.1084
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(lmod.full_faa, ~ . -duration-no_pasg ), test = "Chi")## Single term deletions
##
## Model:
## long.landing ~ aircraft1 + speed_ground + height + pitch
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.58 61.58
## aircraft1 1 86.06 94.06 34.48 4.315e-09 ***
## speed_ground 1 583.66 591.66 532.08 < 2.2e-16 ***
## height 1 75.18 83.18 23.60 1.188e-06 ***
## pitch 1 54.40 62.40 2.82 0.09301 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(lmod.full_faa, ~ . -duration-no_pasg-pitch ), test = "Chi")## Single term deletions
##
## Model:
## long.landing ~ aircraft1 + speed_ground + height
## Df Deviance AIC LRT Pr(>Chi)
## <none> 54.40 62.40
## aircraft1 1 95.06 101.06 40.66 1.814e-10 ***
## speed_ground 1 583.73 589.73 529.33 < 2.2e-16 ***
## height 1 78.16 84.16 23.76 1.090e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Inference: Using drop1 function we get the same model as in above step
Automatic Variable selection Based on AIC
Step 6. Use the R function “Step” to perform forward variable selection using AIC. Compare the result with the table obtained in Step 3. Are the results consistent?
#Stepwise
#create null & full model
lmod.null_faa=lm(long.landing~1, data=FAAlog)
#full model
lmod.full_faa<-glm(long.landing~aircraft1+duration+no_pasg+speed_ground+height+pitch,family=binomial,FAAlog)
model_step_F_AIC_LD <- step(lmod.null_faa, scope=list(lower=lmod.null_faa,
upper=lmod.full_faa),
direction='forward')## Start: AIC=-1716.37
## long.landing ~ 1
##
## Df Sum of Sq RSS AIC
## + speed_ground 1 34.043 53.185 -2101.8
## + aircraft1 1 1.532 85.696 -1728.2
## + pitch 1 0.284 86.944 -1716.9
## <none> 87.229 -1716.4
## + height 1 0.048 87.181 -1714.8
## + duration 1 0.027 87.201 -1714.6
## + no_pasg 1 0.024 87.205 -1714.6
##
## Step: AIC=-2101.76
## long.landing ~ speed_ground
##
## Df Sum of Sq RSS AIC
## + aircraft1 1 2.42022 50.765 -2136.2
## + pitch 1 0.74167 52.444 -2110.8
## + height 1 0.24065 52.945 -2103.3
## <none> 53.185 -2101.8
## + no_pasg 1 0.03285 53.153 -2100.2
## + duration 1 0.01861 53.167 -2100.0
##
## Step: AIC=-2136.23
## long.landing ~ speed_ground + aircraft1
##
## Df Sum of Sq RSS AIC
## + height 1 0.27115 50.494 -2138.4
## <none> 50.765 -2136.2
## + pitch 1 0.10301 50.662 -2135.8
## + duration 1 0.04338 50.722 -2134.9
## + no_pasg 1 0.01837 50.747 -2134.5
##
## Step: AIC=-2138.42
## long.landing ~ speed_ground + aircraft1 + height
##
## Df Sum of Sq RSS AIC
## <none> 50.494 -2138.4
## + pitch 1 0.090011 50.404 -2137.8
## + duration 1 0.041892 50.452 -2137.1
## + no_pasg 1 0.024206 50.470 -2136.8Observation: Model suggested by step function is same as above long.landing ~ speed_ground + aircraft1 + height
Based on BIC
Step 7. Use the R function “Step” to perform forward variable selection using BIC. Compare the result with that from the previous step.
#This model auto uses only 195 rows that have no missing values
model_step_F_BIC_LD <- step(lmod.null_faa, scope=list(lower=lmod.null_faa,
upper=lmod.full_faa),
direction='forward',
criterion = "BIC",
k=log(nrow(FAAlog))
)## Start: AIC=-1711.71
## long.landing ~ 1
##
## Df Sum of Sq RSS AIC
## + speed_ground 1 34.043 53.185 -2092.4
## + aircraft1 1 1.532 85.696 -1718.9
## <none> 87.229 -1711.7
## + pitch 1 0.284 86.944 -1707.6
## + height 1 0.048 87.181 -1705.5
## + duration 1 0.027 87.201 -1705.3
## + no_pasg 1 0.024 87.205 -1705.3
##
## Step: AIC=-2092.43
## long.landing ~ speed_ground
##
## Df Sum of Sq RSS AIC
## + aircraft1 1 2.42022 50.765 -2122.2
## + pitch 1 0.74167 52.444 -2096.8
## <none> 53.185 -2092.4
## + height 1 0.24065 52.945 -2089.3
## + no_pasg 1 0.03285 53.153 -2086.2
## + duration 1 0.01861 53.167 -2086.0
##
## Step: AIC=-2122.24
## long.landing ~ speed_ground + aircraft1
##
## Df Sum of Sq RSS AIC
## <none> 50.765 -2122.2
## + height 1 0.27115 50.494 -2119.8
## + pitch 1 0.10301 50.662 -2117.2
## + duration 1 0.04338 50.722 -2116.2
## + no_pasg 1 0.01837 50.747 -2115.9Observation:
Model suggested by step function BIC is not same as above It only selects 2 variables
long.landing ~ speed_ground + aircraft1
HEIGHT : BIC being parsimonious in its approach does not select height as a significant predictor, We can use aliklihood ratio test to find out if adding this variable can improve the overall model
Liklihood ratio test
#Selected model
lmod.s_ahsg<-glm(long.landing~aircraft1+height+speed_ground,family=binomial,FAAlog)
summary(lmod.s_ahsg)##
## Call:
## glm(formula = long.landing ~ aircraft1 + height + speed_ground,
## family = binomial, data = FAAlog)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.35721 -0.00160 -0.00001 0.00000 2.55053
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -98.86483 18.64743 -5.302 1.15e-07 ***
## aircraft1 4.91354 1.10438 4.449 8.62e-06 ***
## height 0.22063 0.06028 3.660 0.000252 ***
## speed_ground 0.88939 0.16684 5.331 9.79e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 598.240 on 782 degrees of freedom
## Residual deviance: 54.401 on 779 degrees of freedom
## AIC: 62.401
##
## Number of Fisher Scoring iterations: 11#BIC Model
lmod.asg<-glm(long.landing~aircraft1+speed_ground,family=binomial,FAAlog)
summary(lmod.asg)##
## Call:
## glm(formula = long.landing ~ aircraft1 + speed_ground, family = binomial,
## data = FAAlog)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.27583 -0.01388 -0.00038 0.00000 2.58037
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -61.01550 8.91176 -6.847 7.56e-12 ***
## aircraft1 3.27823 0.74478 4.402 1.07e-05 ***
## speed_ground 0.58743 0.08668 6.777 1.23e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 598.240 on 782 degrees of freedom
## Residual deviance: 78.164 on 780 degrees of freedom
## AIC: 84.164
##
## Number of Fisher Scoring iterations: 10anova(lmod.asg,lmod.s_ahsg,test="Chi")## Analysis of Deviance Table
##
## Model 1: long.landing ~ aircraft1 + speed_ground
## Model 2: long.landing ~ aircraft1 + height + speed_ground
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 780 78.164
## 2 779 54.401 1 23.762 1.09e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Inference:
Test suggests that bigger model model is significantly different from smaller model hence we select the bigger model.
Final Analysis: Report
Step 8. Now,You are scheduled to meet with an FAA agent who wants to know “what are risk factors for long landings and how do they influence its occurrence?”.
#AIC BIC Comparision
print(paste("AIC of BIC suggested (aircraft+speed_air) model is",round(AIC(lmod.asg),4)),quote=FALSE)## [1] AIC of BIC suggested (aircraft+speed_air) model is 84.1637print(paste("BIC of BIC suggested (aircraft+speed_air) model is",round(BIC(lmod.asg),4)),quote=FALSE)## [1] BIC of BIC suggested (aircraft+speed_air) model is 98.1531print(paste("AIC of selected model(aircraft+speed_ground+height) is",round(AIC(lmod.s_ahsg),4)),quote=FALSE)## [1] AIC of selected model(aircraft+speed_ground+height) is 62.4014print(paste("BIC of selected model(aircraft+speed_ground+height) is",round(BIC(lmod.s_ahsg),4)),quote=FALSE)## [1] BIC of selected model(aircraft+speed_ground+height) is 81.0539f<-cbind(coef(summary(lmod.s_ahsg)),exp(coef(summary(lmod.s_ahsg))[,1]))
colnames(f)<-c("Estimate","Std. Error","z value","Pr(>|z|)","Odds Ratio" )
f## Estimate Std. Error z value Pr(>|z|) Odds Ratio
## (Intercept) -98.8648279 18.6474315 -5.301793 1.146706e-07 1.157579e-43
## aircraft1 4.9135434 1.1043832 4.449129 8.621918e-06 1.361209e+02
## height 0.2206339 0.0602816 3.660054 2.521616e-04 1.246867e+00
## speed_ground 0.8893937 0.1668449 5.330661 9.785594e-08 2.433654e+00We observe that
Sample size: We started with 783 line items in dataset, and we have been able to maintain the same sample size for the analysis of long landing Relevent predictors: aircraft ,speed_ground height Missing values: Only the speed_air column has 583(75%) missing values, we have not used it in out analysis.
Inferences: * Three variables aircraft, speed_air, height significantly effect the long.landing binary Response variable * The AIC, BIC (Criterian of model selection) are lowest when we include only these 3 predictors. * Estimate of aircraft variable is highest 4.9 followed by speed ground : .9 and then the lowest is of height: 0.22 * As the predictor increase by 1 unit the log odds increase linearly by estimate coefficient value. This also means that odds of success(long.landing to happen) also increases (not linearly but does increase) exp(coefficient value) value * Binary Predictor: if its a boeing plane the log odds are higher by a value of 4.9, because log odds are propotional to odds(although not lineraly), hence odds or risky landing are higher with a boeing plane.
Hence these predictors are very important to correctly predic the landing distance over run.
4.2 Risky Landing
Identifying important factors using the binary data of “risky.landing”
Step 9. Repeat Steps 1-7 but using “risky.landing” as the binary response
Univariate Analysis
Step2.b
Identifying important factors using the binary data of “risky.landing”.
Use a pie chart or a histogram to show the distribution of “risky.landing”.
FAAlog$aircraft1<-ifelse(FAAlog$aircraft=='boeing',as.numeric(1),as.numeric(0))
par(mfrow=c(1,2))
plot(risky.landing~speed_ground,FAAlog)
plot(jitter(risky.landing,0.2)~jitter(speed_ground),FAAlog,xlab ="speed ground",ylab ="long.landing",pch=".")
par(mfrow=c(1,2))
plot(risky.landing~speed_air,FAAlog)
plot(jitter(risky.landing,0.2)~jitter(speed_air),FAAlog,xlab ="speed_air",ylab ="long.landing",pch=".")
par(mfrow=c(1,2))
plot(long.landing~aircraft1,FAAlog)
plot(jitter(long.landing,0.2)~jitter(aircraft1),FAAlog,xlab ="aircraft",ylab ="risky.landing",pch=".")
#names(FAA_U_RmvNA)
#FAAlog_R<-FAA_U_RmvNA[,c(1,2,3,4,5,6,7,10)]
# A model with all variables available
R_lmod.aircraft<-glm(risky.landing~aircraft,family=binomial,FAAlog)
R_lmod.duration<-glm(risky.landing~duration,family=binomial,FAAlog)
R_lmod.no_pasg<-glm(risky.landing~no_pasg,family=binomial,FAAlog)
R_lmod.speed_ground<-glm(risky.landing~speed_ground,family=binomial,FAAlog)
R_lmod.speed_air<-glm(risky.landing~speed_air,family=binomial,FAAlog)
R_lmod.height<-glm(risky.landing~height,family=binomial,FAAlog)
R_lmod.pitch<-glm(risky.landing~pitch,family=binomial,FAAlog)
Rv1<-c('aircraft',round(coef(R_lmod.aircraft)[1],2),coef(R_lmod.aircraft)[2],exp(coef(R_lmod.aircraft))[1],exp(coef(R_lmod.aircraft))[2],
sign<-ifelse(coef(R_lmod.aircraft)[2]==-1,"Negative","Positive"),coef(summary(R_lmod.aircraft))[,4])
Rv2<-c('duration',coef(R_lmod.duration)[1],coef(R_lmod.duration)[2],exp(coef(R_lmod.duration))[1],exp(coef(R_lmod.duration))[2],
sign<-ifelse(coef(R_lmod.duration)[2]==-1,"Negative","Positive"),coef(summary(R_lmod.duration))[,4])
Rv3<-c('no_pasg',coef(R_lmod.no_pasg)[1],coef(R_lmod.no_pasg)[2],exp(coef(R_lmod.no_pasg))[1],exp(coef(R_lmod.no_pasg))[2],
sign<-ifelse(sign(coef(R_lmod.no_pasg)[2])==-1,"Negative","Positive"),coef(summary(R_lmod.no_pasg))[,4])
Rv4<-c('speed_ground',coef(R_lmod.speed_ground)[1],coef(R_lmod.speed_ground)[2],exp(coef(R_lmod.speed_ground))[1],exp(coef(R_lmod.speed_ground))[2],
sign<-ifelse(coef(R_lmod.speed_ground)[2]==-1,"Negative","Positive"),coef(summary(R_lmod.speed_ground))[,4])
Rv5<-c('speed_air',coef(R_lmod.speed_air)[1],coef(R_lmod.speed_air)[2],exp(coef(R_lmod.speed_air))[1],exp(coef(R_lmod.speed_air))[2],
sign<-ifelse(coef(R_lmod.speed_air)[2]==-1,"Negative","Positive"),coef(summary(R_lmod.speed_air))[,4])
Rv6<-c('height',coef(R_lmod.height)[1],coef(R_lmod.height)[2],exp(coef(R_lmod.height))[1],exp(coef(R_lmod.height))[2],
sign<-ifelse(coef(R_lmod.height)[2]==-1,"Negative","Positive"),coef(summary(R_lmod.height))[,4])
Rv7<-c('pitch',coef(R_lmod.pitch)[1],coef(R_lmod.pitch)[2],exp(coef(R_lmod.pitch))[1],exp(coef(R_lmod.pitch))[2],
sign<-ifelse(coef(R_lmod.pitch)[2]==-1,"Negative","Positive"),coef(summary(R_lmod.pitch))[,4])
Rd<-data.frame(rbind(Rv1,Rv2,Rv3,Rv4,Rv5,Rv6,Rv7))
colnames(Rd)<-c('Name','B0','B1','OddsRto_B0','OddsRto_B1','Sign','Pvalue_B0','Pvalue_B1')
Rd$B0<-round(as.numeric(as.character(Rd$B0)),2)
Rd$B1<-round(as.numeric(as.character(Rd$B1)),2)
Rd$OddsRto_B0<-round(as.numeric(as.character(Rd$OddsRto_B0)),2)
Rd$OddsRto_B1<-round(as.numeric(as.character(Rd$OddsRto_B1)),2)
Rd$Pvalue_B0<-as.numeric(as.character(Rd$Pvalue_B0))
Rd$Pvalue_B1<-as.numeric(as.character(Rd$Pvalue_B1))
Rd<-arrange(Rd,Pvalue_B1,Pvalue_B0)
Rd$signif<-ifelse(Rd$Pvalue_B1<0.05,'S','I')
Rd[,c(1,2,3,6,7,8,9)]## Name B0 B1 Sign Pvalue_B0 Pvalue_B1 signif
## 1 speed_ground -64.40 0.60 Positive 9.638076e-08 1.085926e-07 S
## 2 speed_air -91.48 0.85 Positive 5.733858e-06 5.720457e-06 S
## 3 aircraft -3.04 0.93 Positive 2.175652e-36 1.455107e-03 S
## 4 no_pasg -1.00 -0.03 Negative 3.476978e-01 1.594340e-01 I
## 5 pitch -3.83 0.33 Positive 3.336993e-04 2.018700e-01 I
## 6 duration -2.31 0.00 Positive 2.428451e-07 6.730755e-01 I
## 7 height -2.37 0.00 Positive 6.269855e-08 7.777416e-01 IObservation:
- 3 predictors seem to be significant based on pvalue analysis, namely: speed_ground,speed_air,aircraft
- All coefficients are positive: which means they positively effect the landing distance
- Aircraft has the biggest coefficient follwed by speed air and speed ground.
- We also observe that pvalue for height is not significant.
** Fit a model**
Here we fit a full model, based on prior collinearity analysis we remove speed air fit a full model after removing collinerity
#full model
Rlmod.full_faa<-glm(risky.landing~aircraft1+duration+no_pasg+speed_ground+height+pitch,family=binomial,FAAlog)
summary(Rlmod.full_faa)#aircraft and speed_air##
## Call:
## glm(formula = risky.landing ~ aircraft1 + duration + no_pasg +
## speed_ground + height + pitch, family = binomial, data = FAAlog)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.44763 -0.00011 0.00000 0.00000 1.85688
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.022e+02 2.811e+01 -3.635 0.000278 ***
## aircraft1 4.406e+00 1.562e+00 2.821 0.004783 **
## duration 7.386e-04 1.214e-02 0.061 0.951498
## no_pasg -8.590e-02 6.011e-02 -1.429 0.152955
## speed_ground 9.366e-01 2.476e-01 3.782 0.000155 ***
## height 4.214e-02 4.618e-02 0.913 0.361502
## pitch 6.083e-01 8.000e-01 0.760 0.447088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 423.535 on 782 degrees of freedom
## Residual deviance: 36.372 on 776 degrees of freedom
## AIC: 50.372
##
## Number of Fisher Scoring iterations: 12drop1 analysis:
drop1(Rlmod.full_faa,test ="Chi")## Single term deletions
##
## Model:
## risky.landing ~ aircraft1 + duration + no_pasg + speed_ground +
## height + pitch
## Df Deviance AIC LRT Pr(>Chi)
## <none> 36.37 50.37
## aircraft1 1 49.90 61.90 13.53 0.0002345 ***
## duration 1 36.38 48.38 0.00 0.9515067
## no_pasg 1 38.62 50.62 2.25 0.1334514
## speed_ground 1 410.63 422.63 374.26 < 2.2e-16 ***
## height 1 37.23 49.23 0.86 0.3549199
## pitch 1 36.96 48.96 0.59 0.4423665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(Rlmod.full_faa, ~ . -duration), test = "Chi")## Single term deletions
##
## Model:
## risky.landing ~ aircraft1 + no_pasg + speed_ground + height +
## pitch
## Df Deviance AIC LRT Pr(>Chi)
## <none> 36.38 48.38
## aircraft1 1 49.91 59.91 13.53 0.0002347 ***
## no_pasg 1 38.82 48.82 2.44 0.1182167
## speed_ground 1 410.75 420.75 374.37 < 2.2e-16 ***
## height 1 37.30 47.30 0.93 0.3351600
## pitch 1 36.98 46.98 0.61 0.4352093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(Rlmod.full_faa, ~ . -duration-pitch ), test = "Chi")## Single term deletions
##
## Model:
## risky.landing ~ aircraft1 + no_pasg + speed_ground + height
## Df Deviance AIC LRT Pr(>Chi)
## <none> 36.98 46.98
## aircraft1 1 57.06 65.06 20.08 7.439e-06 ***
## no_pasg 1 39.30 47.30 2.31 0.1285
## speed_ground 1 410.77 418.77 373.78 < 2.2e-16 ***
## height 1 37.56 45.56 0.57 0.4485
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(Rlmod.full_faa, ~ . -duration-pitch-height ), test = "Chi")## Single term deletions
##
## Model:
## risky.landing ~ aircraft1 + no_pasg + speed_ground
## Df Deviance AIC LRT Pr(>Chi)
## <none> 37.56 45.56
## aircraft1 1 57.18 63.18 19.62 9.45e-06 ***
## no_pasg 1 39.96 45.96 2.40 0.1216
## speed_ground 1 410.79 416.79 373.24 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(update(Rlmod.full_faa, ~ . -duration-height-pitch-no_pasg), test = "Chi")## Single term deletions
##
## Model:
## risky.landing ~ aircraft1 + speed_ground
## Df Deviance AIC LRT Pr(>Chi)
## <none> 39.96 45.96
## aircraft1 1 57.99 61.99 18.03 2.171e-05 ***
## speed_ground 1 412.53 416.53 372.57 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#suggested model : aircraft and speed ground only
Rlmod.s<-glm(risky.landing~aircraft1+speed_ground,family=binomial,FAAlog)
summary(Rlmod.s)##
## Call:
## glm(formula = risky.landing ~ aircraft1 + speed_ground, family = binomial,
## data = FAAlog)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.22465 -0.00014 0.00000 0.00000 1.60326
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -100.6448 24.8834 -4.045 5.24e-05 ***
## aircraft1 3.9763 1.2520 3.176 0.00149 **
## speed_ground 0.9132 0.2258 4.044 5.25e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 423.535 on 782 degrees of freedom
## Residual deviance: 39.955 on 780 degrees of freedom
## AIC: 45.955
##
## Number of Fisher Scoring iterations: 12Automatic Variable Selection
AIC Step 6(2). Use the R function “Step” to perform forward variable selection using AIC. Compare the result with the table obtained in Step 3.
Are the results consistent?
#Stepwise
#create null & full model
Rlmod.null_faa=lm(risky.landing~1, data=FAAlog)
#full model
Rlmod.full_faa<-glm(risky.landing~aircraft1+duration+no_pasg+speed_ground+height+pitch,family=binomial,FAAlog)
#This model auto uses only 195 rows that have no missing values
Rmodel_step_F_AIC <- step(Rlmod.null_faa, scope=list(lower=Rlmod.null_faa,
upper=Rlmod.full_faa),
direction='forward')## Start: AIC=-2071.78
## risky.landing ~ 1
##
## Df Sum of Sq RSS AIC
## + speed_ground 1 16.5163 38.886 -2347.0
## + aircraft1 1 0.7593 54.643 -2080.6
## <none> 55.402 -2071.8
## + no_pasg 1 0.1403 55.262 -2071.8
## + pitch 1 0.1155 55.287 -2071.4
## + duration 1 0.0126 55.390 -2070.0
## + height 1 0.0056 55.397 -2069.9
##
## Step: AIC=-2346.96
## risky.landing ~ speed_ground
##
## Df Sum of Sq RSS AIC
## + aircraft1 1 1.19415 37.692 -2369.4
## + pitch 1 0.32298 38.563 -2351.5
## + no_pasg 1 0.15469 38.731 -2348.1
## <none> 38.886 -2347.0
## + height 1 0.01303 38.873 -2345.2
## + duration 1 0.00953 38.876 -2345.2
##
## Step: AIC=-2369.38
## risky.landing ~ speed_ground + aircraft1
##
## Df Sum of Sq RSS AIC
## + no_pasg 1 0.130547 37.561 -2370.1
## <none> 37.692 -2369.4
## + pitch 1 0.034664 37.657 -2368.1
## + duration 1 0.021931 37.670 -2367.8
## + height 1 0.018309 37.673 -2367.8
##
## Step: AIC=-2370.1
## risky.landing ~ speed_ground + aircraft1 + no_pasg
##
## Df Sum of Sq RSS AIC
## <none> 37.561 -2370.1
## + pitch 1 0.033941 37.527 -2368.8
## + height 1 0.022276 37.539 -2368.6
## + duration 1 0.018207 37.543 -2368.5#aic model
Rlmod.AIC<-glm(risky.landing ~ speed_ground + aircraft1 + no_pasg,family=binomial,FAAlog)AIC STEP function suggests that significant predictors are:
risky.landing ~ speed_ground + aircraft1 + no_pasg
BIC
Step 7(2). Use the R function “Step” to perform forward variable selection using BIC. Compare the result with that from the previous step.
#This model auto uses only 195 rows that have no missing values
Rmodel_step_F_BIC <- step(Rlmod.null_faa, scope=list(lower=Rlmod.null_faa,
upper=Rlmod.full_faa),
direction='forward',
criterion = "BIC",
k=log(nrow(FAAlog))
)## Start: AIC=-2067.12
## risky.landing ~ 1
##
## Df Sum of Sq RSS AIC
## + speed_ground 1 16.5163 38.886 -2337.6
## + aircraft1 1 0.7593 54.643 -2071.3
## <none> 55.402 -2067.1
## + no_pasg 1 0.1403 55.262 -2062.4
## + pitch 1 0.1155 55.287 -2062.1
## + duration 1 0.0126 55.390 -2060.6
## + height 1 0.0056 55.397 -2060.5
##
## Step: AIC=-2337.63
## risky.landing ~ speed_ground
##
## Df Sum of Sq RSS AIC
## + aircraft1 1 1.19415 37.692 -2355.4
## <none> 38.886 -2337.6
## + pitch 1 0.32298 38.563 -2337.5
## + no_pasg 1 0.15469 38.731 -2334.1
## + height 1 0.01303 38.873 -2331.2
## + duration 1 0.00953 38.876 -2331.2
##
## Step: AIC=-2355.39
## risky.landing ~ speed_ground + aircraft1
##
## Df Sum of Sq RSS AIC
## <none> 37.692 -2355.4
## + no_pasg 1 0.130547 37.561 -2351.4
## + pitch 1 0.034664 37.657 -2349.4
## + duration 1 0.021931 37.670 -2349.2
## + height 1 0.018309 37.673 -2349.1Model suggested by BIC is: risky.landing ~ speed_ground + aircraft1
Liklihood Ratio test
#model
anova(Rlmod.s,Rlmod.AIC,test="Chi")#cannot reject null hypothesis## Analysis of Deviance Table
##
## Model 1: risky.landing ~ aircraft1 + speed_ground
## Model 2: risky.landing ~ speed_ground + aircraft1 + no_pasg
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 780 39.955
## 2 779 37.559 1 2.3964 0.1216Interpretation
Liklihood ratio test suggests we can select the smaller model Hence, we select the model with Speed_ground and aircraft
#AIC MODEL
print(paste("AIC of AIC model is",round(AIC(Rlmod.AIC),4)),quote=FALSE)## [1] AIC of AIC model is 45.5589print(paste("AIC of AIC model is",round(BIC(Rlmod.AIC),4)),quote=FALSE)## [1] AIC of AIC model is 64.2115#SELECTED FINAL BIC & DROP FUNCTION
print(paste("AIC of FINAL model(DROP FUNCTION) is",round(AIC(Rlmod.s),4)),quote=FALSE)## [1] AIC of FINAL model(DROP FUNCTION) is 45.9553print(paste("AIC of FINAL model(DROP FUNCTION) is",round(BIC(Rlmod.s),4)),quote=FALSE)## [1] AIC of FINAL model(DROP FUNCTION) is 59.9447coef(summary(Rlmod.s))## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -100.6447977 24.8834111 -4.044654 5.240037e-05
## aircraft1 3.9763497 1.2520073 3.175980 1.493315e-03
## speed_ground 0.9131614 0.2258071 4.043989 5.254930e-05Rf<-cbind(coef(summary(Rlmod.s)),exp(coef(summary(Rlmod.s))[,1]))
colnames(Rf)<-c("Estimate","Std. Error","z value","Pr(>|z|)","Odds Ratio" )
Rf## Estimate Std. Error z value Pr(>|z|) Odds Ratio
## (Intercept) -100.6447977 24.8834111 -4.044654 5.240037e-05 1.952179e-44
## aircraft1 3.9763497 1.2520073 3.175980 1.493315e-03 5.332204e+01
## speed_ground 0.9131614 0.2258071 4.043989 5.254930e-05 2.492189e+00We observe that Sample size: We started with 783 line items in dataset, and we have been able to maintain the same sample size for the analysis of long landing Relevent predictors: aircraft ,speed_ground Missing values: Only the speed_air column has 583(75%) missing values, we have not used it in out analysis.
Inferences: * Two variables aircraft, speed_air significantly effect the risky.landing binary Response variable * The AIC is marginally higher but BIC is loweres for this model * Estimate of aircraft variable is highest 4.0(positive) followed by speed ground .91(positive). * Continious predictor: As the predictor increase by 1 unit the log odds increase linearly by estimate coefficient value which inturn means odds of success(risky landing to happen) also increases (not linearly but does increase)exp(coefficient value). * Binary Predictor: if its a boeing plane the log odds are higher by a value of 4.0, because log odds are propotional to odds(although not lineraly), hence odds or risky landing are higher with a boeing plane.
Hence these predictors are very important to correctly predic the landing distance over run.
4.3 Model Comparision
Step 11. Compare the two models built for “long.landing” and “risky.landing”
Model selection long landing model uses height, aircraft and speed_air predictors where as risky landing uses only aircraft and speed_air as predictors. For predicting long landing height is and important predictor but its not for risky landing. Apparently if Landing distance increases more that 3000 (risky landing), height does not play an important role in predicting risky landing. Coefficient of Aircraft: long landing Beta0 =4.9 and risky landing Beta0=4.0 For this case Boeing = 1 because log odds increases linearly with Beta0, which inturn in propotional to Odds. Hence,if boeing planes are selected then there higher odds for long landing than risky landing.
4.4 ROC Curve
Step 12. Plot the ROC curve (sensitivity versus 1-specificity) for each model (see pp.32-33 in Lecture 4 slides). Draw the two curves in the same plot. Do you have any comment?
long landing
linpred.l<-predict(lmod.s_ahsg) # Linear predictor
predprob.l<-predict(lmod.s_ahsg,type ="response") #predicted probabilities
predout.l<-(ifelse(predprob.l<0.5,"no","yes")) ### Predicted outcomes using 0.5 as the threshold
FAAlog.l<-data.frame(FAAlog,predprob.l,predout.l)
xtabs(~long.landing+predout.l,FAAlog.l)## predout.l
## long.landing no yes
## 0 679 4
## 1 6 94thresh<-seq(0.01,0.5,0.01)
sensitivity.l<-specificity.l<-rep(NA,length(thresh))
for(j in seq (along=thresh)){
pp.l<-ifelse(FAAlog.l$predprob.l< thresh[j],"no","yes")
xx.l<-xtabs(~long.landing+pp.l,FAAlog.l)
specificity.l[j]<-xx.l[1,1]/(xx.l[1,1]+xx.l[1,2])
sensitivity.l[j]<-xx.l[2,2]/(xx.l[2,1]+xx.l[2,2])
}
par(mfrow=c(1,2))
matplot(thresh,cbind(sensitivity.l,specificity.l),type="l", xlab ="Threshold",ylab="proportion",lty=1:2)
plot(1-specificity.l,sensitivity.l,type="l",col="Blue");abline(0,1,lty=2)
Risky landing
linpredF<-predict(Rlmod.s) # Linear predictor
predprobF<-predict(Rlmod.s,type ="response") #predicted probabilities
predoutF<-(ifelse(predprobF<0.5,"no","yes")) ### Predicted outcomes using 0.5 as the threshold
FAAlog.F<-data.frame(FAAlog,predprobF,predoutF)
xtabs(~risky.landing+predoutF,FAAlog)## predoutF
## risky.landing no yes
## 0 717 6
## 1 6 54threshF<-seq(0.01,0.999,0.01)
sensitivityF<-specificityF<-rep(NA,length(threshF))
for(j in seq (along=threshF)){
ppF<-ifelse(FAAlog.F$predprobF< threshF[j],"no","yes")
xxF<-xtabs(~risky.landing+ppF,FAAlog)
specificityF[j]<-xxF[1,1]/(xxF[1,1]+xxF[1,2])
sensitivityF[j]<-xxF[2,2]/(xxF[2,1]+xxF[2,2])
}
par(mfrow=c(1,2))
matplot(threshF,cbind(sensitivityF,specificityF),type="l", xlab ="Threshold",ylab="proportion",lty=1:2)
plot(1-specificityF,sensitivityF,type="l",col="Red");abline(0,1,lty=2)
PLOTING both togather:
plot(1-specificityF,sensitivityF,type="l",col="Red",xlab="1-Specificity",ylab="Sensitivity")
lines(1-specificity.l,sensitivity.l,col="Blue")
4.5 Prediction
Step 13.
A commercial airplane is passing over the threshold of the runway, at this moment we have its basic information and measures of its airborne performance (Boeing, duration=200, no_pasg=80, speed_ground=115, speed_air=120, height=40, pitch=4). Predict its probability of being a long landing and a risky landing, respectively. Report the predicted probability as well as its 95% confidence interval. Long landing
new.ind<-data.frame(aircraft1=1,duration=200,no_pasg=80,speed_ground=115,speed_air=120,height=40,pitch=4)
predict(lmod.s_ahsg,newdata=new.ind,type ="link") ### Linear predictor## 1
## 17.15434predict(lmod.s_ahsg,newdata= new.ind,type ="response") ### Probability## 1
## 1### Predict the value with its standard error
predict(lmod.s_ahsg,newdata=new.ind,type ="link",se=T) ### Linear predictor## $fit
## 1
## 17.15434
##
## $se.fit
## [1] 3.366708
##
## $residual.scale
## [1] 1predict(lmod.s_ahsg,newdata= new.ind,type ="response",se=T) ### Probability## $fit
## 1
## 1
##
## $se.fit
## 1
## 1.194452e-07
##
## $residual.scale
## [1] 1round(ilogit(c( 17.15434+1.96*3.366708, 17.15434+1.96*3.366708)),3 ) ### Confidence interval## [1] 1 1round(c(1+1.96*1.194452e-07,1-1.96*1.194452e-07), 3 ) ### Confidence interval## [1] 1 1We observe that confidence Interval is 1,1 which is an indicator of better prediction ability of our selected model
Risky landing
predict(Rlmod.s,newdata=new.ind,type ="link") ### Linear predictor## 1
## 8.345109predict(Rlmod.s,newdata= new.ind,type ="response") ### Probability## 1
## 0.9997625### Predict the value with its standard error
predict(Rlmod.s,newdata=new.ind,type ="link",se=T) ### Linear predictor## $fit
## 1
## 8.345109
##
## $se.fit
## [1] 2.098615
##
## $residual.scale
## [1] 1predict(Rlmod.s,newdata= new.ind,type ="response",se=T) ### Probability## $fit
## 1
## 0.9997625
##
## $se.fit
## 1
## 0.0004983008
##
## $residual.scale
## [1] 1round(ilogit(c( 8.345109+1.96*2.098615, 8.345109+1.96*2.098615)),5 ) ### Confidence interval## [1] 1 1round(c(0.9997625+1.96*0.0004983008 ,0.9997625-1.96*0.0004983008), 5 ) ### Confidence interval## [1] 1.00074 0.99879We observe that confidence Interval marginally different ,because the difference is very small we can avoid it.
5.0 Other functions
5.1 Probit vz Cloglog
Compare models with different link functions
Step 14. For the binary response “risky landing”, fit the following models using the risk factors identified in Steps 9-10: Probit model Hazard model with complementary log-log link Compare these two models with the logistic model. Do you have any comments?
lmod.s_ahsg_probit<-glm(long.landing~aircraft1+height+speed_ground,family=binomial(link=probit),FAAlog)
lmod.s_ahsg_cloglog<-glm(long.landing~aircraft1+height+speed_ground,family=binomial(link=cloglog),FAAlog)
df.compare.coef<-cbind(round(coef(lmod.s_ahsg),3),round(coef(lmod.s_ahsg_probit),3),round(coef(lmod.s_ahsg_cloglog),3))
colnames(df.compare.coef)<-c('logit','probit','cloglog')
df.compare.coef## logit probit cloglog
## (Intercept) -98.865 -53.932 -67.868
## aircraft1 4.914 2.694 3.664
## height 0.221 0.117 0.153
## speed_ground 0.889 0.486 0.602Inference: coefficient of variable is highest in logit function followed by comp log log function and then by probit
predval<-sapply(list(lmod.s_ahsg,lmod.s_ahsg_probit,lmod.s_ahsg_cloglog),fitted)
colnames(predval)<- c('logit','probit','cloglog')Compare top ten values
predval[1:10 ,]## logit probit cloglog
## 1 9.996897e-01 9.999947e-01 1.000000e+00
## 2 9.305530e-01 9.202875e-01 9.652767e-01
## 3 2.649443e-12 2.220446e-16 8.265371e-09
## 4 1.947901e-05 1.395946e-09 3.829680e-04
## 5 2.220446e-16 2.220446e-16 8.322862e-11
## 6 1.323259e-08 2.220446e-16 2.878625e-06
## 7 2.220446e-16 2.220446e-16 8.654803e-13
## 8 2.220446e-16 2.220446e-16 2.117752e-12
## 9 3.894756e-05 1.207688e-08 6.235964e-04
## 10 2.220446e-16 2.220446e-16 5.111611e-10Compare values between 03 and 0.5
predval[fitted(lmod.s_ahsg)>0.3 & fitted(lmod.s_ahsg) <0.5,]## logit probit cloglog
## 24 0.3968172 0.4050834 0.3547947
## 60 0.3434765 0.3422366 0.3271934
## 202 0.4178148 0.4349099 0.3606496
## 399 0.3284088 0.3449841 0.2215693
## 667 0.4584596 0.4561809 0.3081406
## 693 0.3688203 0.3523023 0.2688676Compare Values <0.01
predval[fitted(lmod.s_ahsg)<0.01,]## logit probit cloglog
## 3 2.649443e-12 2.220446e-16 8.265371e-09
## 4 1.947901e-05 1.395946e-09 3.829680e-04
## 5 2.220446e-16 2.220446e-16 8.322862e-11
## 6 1.323259e-08 2.220446e-16 2.878625e-06
## 7 2.220446e-16 2.220446e-16 8.654803e-13
## 8 2.220446e-16 2.220446e-16 2.117752e-12
## 9 3.894756e-05 1.207688e-08 6.235964e-04
## 10 2.220446e-16 2.220446e-16 5.111611e-10
## 11 2.220446e-16 2.220446e-16 1.792280e-11
## 14 2.220446e-16 2.220446e-16 2.986759e-10
## 15 2.220446e-16 2.220446e-16 8.686422e-14
## 16 2.171105e-06 4.695932e-13 8.653623e-05
## 17 1.475757e-05 6.057585e-10 3.106324e-04
## 20 2.477480e-10 2.220446e-16 1.858413e-07
## 21 8.130356e-08 2.220446e-16 9.361882e-06
## 23 2.220446e-16 2.220446e-16 1.133932e-10
## 25 6.930428e-10 2.220446e-16 3.835286e-07
## 26 2.393752e-03 4.694756e-04 9.808371e-03
## 27 4.437937e-11 2.220446e-16 5.738253e-08
## 28 1.361982e-11 2.220446e-16 2.784306e-08
## 31 2.220446e-16 2.220446e-16 1.403836e-15
## 32 1.923537e-08 2.220446e-16 3.452152e-06
## 33 2.118848e-09 2.220446e-16 8.395121e-07
## 34 9.772867e-10 2.220446e-16 5.135136e-07
## 36 2.220446e-16 2.220446e-16 3.636170e-11
## 37 2.220446e-16 2.220446e-16 7.045458e-16
## 38 2.499487e-10 2.220446e-16 1.918001e-07
## 39 1.372544e-03 1.711095e-04 6.463745e-03
## 40 2.220446e-16 2.220446e-16 7.702137e-10
## 43 2.220446e-16 2.220446e-16 1.408401e-11
## 44 4.977745e-10 2.220446e-16 2.872163e-07
## 45 1.321884e-12 2.220446e-16 5.364810e-09
## 47 1.001714e-12 2.220446e-16 4.266222e-09
## 48 1.774328e-12 2.220446e-16 6.338201e-09
## 49 7.358823e-06 3.403981e-11 2.118102e-04
## 50 2.749394e-13 2.220446e-16 1.981301e-09
## 51 2.220446e-16 2.220446e-16 4.078216e-10
## 52 1.113192e-06 2.084697e-14 5.898479e-05
## 53 2.220446e-16 2.220446e-16 4.994917e-13
## 54 2.916047e-11 2.220446e-16 4.216338e-08
## 55 2.220446e-16 2.220446e-16 1.642461e-11
## 57 1.803762e-09 2.220446e-16 6.700103e-07
## 58 2.220446e-16 2.220446e-16 1.762804e-12
## 59 2.220446e-16 2.220446e-16 1.907846e-10
## 61 1.952284e-04 1.181228e-06 1.906422e-03
## 62 9.977082e-03 5.912945e-03 2.544436e-02
## 63 2.652202e-05 3.428191e-09 4.842837e-04
## 65 6.884368e-08 2.220446e-16 8.130156e-06
## 66 2.220446e-16 2.220446e-16 9.690998e-15
## 67 2.587148e-06 6.926383e-13 1.030360e-04
## 68 3.352990e-03 9.834535e-04 1.185824e-02
## 69 9.356860e-07 1.651812e-14 4.813305e-05
## 70 2.220446e-16 2.220446e-16 1.026549e-10
## 71 2.172919e-11 2.220446e-16 3.866565e-08
## 73 2.220446e-16 2.220446e-16 1.276561e-13
## 74 2.220446e-16 2.220446e-16 6.925692e-14
## 75 3.819242e-03 1.276356e-03 1.283134e-02
## 76 1.090101e-04 3.483870e-07 1.162850e-03
## 77 2.220446e-16 2.220446e-16 8.383112e-14
## 79 2.220446e-16 2.220446e-16 4.352401e-13
## 80 2.818538e-05 6.345523e-09 4.592297e-04
## 82 4.903992e-09 2.220446e-16 1.408234e-06
## 83 2.220446e-16 2.220446e-16 3.114431e-14
## 84 2.220446e-16 2.220446e-16 1.571936e-11
## 85 2.738599e-10 2.220446e-16 2.065941e-07
## 86 4.095813e-09 2.220446e-16 1.307216e-06
## 88 2.220446e-16 2.220446e-16 1.696709e-14
## 89 2.220446e-16 2.220446e-16 2.220446e-16
## 90 2.954052e-08 2.220446e-16 4.712627e-06
## 92 6.443333e-08 2.220446e-16 9.000060e-06
## 94 1.076996e-09 2.220446e-16 5.243254e-07
## 96 1.404828e-13 2.220446e-16 1.083664e-09
## 98 6.507477e-13 2.220446e-16 3.083018e-09
## 99 2.299161e-07 2.220446e-16 1.878130e-05
## 100 2.220446e-16 2.220446e-16 1.462362e-14
## 101 3.425497e-11 2.220446e-16 5.182286e-08
## 103 2.220446e-16 2.220446e-16 1.117201e-12
## 104 1.515498e-04 8.739729e-07 1.449668e-03
## 105 2.220446e-16 2.220446e-16 2.220446e-16
## 106 7.338978e-08 2.220446e-16 8.346094e-06
## 107 2.375758e-07 2.220446e-16 1.996669e-05
## 108 2.668875e-07 2.220446e-16 1.972182e-05
## 109 1.803074e-05 9.513120e-10 3.734113e-04
## 110 2.220446e-16 2.220446e-16 1.535181e-11
## 111 2.220446e-16 2.220446e-16 5.943733e-12
## 112 1.058865e-13 2.220446e-16 9.753198e-10
## 113 1.208937e-11 2.220446e-16 2.407592e-08
## 116 4.831513e-03 1.784207e-03 1.555917e-02
## 117 2.220446e-16 2.220446e-16 6.734370e-10
## 118 3.971879e-13 2.220446e-16 2.514009e-09
## 121 2.417093e-03 4.180791e-04 1.037614e-02
## 122 7.345929e-06 3.273395e-11 2.129407e-04
## 123 2.220446e-16 2.220446e-16 2.366008e-10
## 125 2.220446e-16 2.220446e-16 3.973961e-15
## 126 2.220446e-16 2.220446e-16 4.816871e-12
## 127 2.331133e-12 2.220446e-16 7.773739e-09
## 128 1.587999e-06 1.554755e-13 6.816645e-05
## 129 1.272344e-10 2.220446e-16 1.200278e-07
## 130 8.531027e-12 2.220446e-16 2.059146e-08
## 131 2.220446e-16 2.220446e-16 8.140210e-13
## 132 2.220446e-16 2.220446e-16 5.357042e-11
## 133 2.453217e-05 2.890499e-09 4.509285e-04
## 134 1.295993e-10 2.220446e-16 1.119064e-07
## 135 2.220446e-16 2.220446e-16 5.544062e-13
## 137 3.891373e-05 1.244378e-08 6.185466e-04
## 138 1.882958e-11 2.220446e-16 3.147879e-08
## 139 2.220446e-16 2.220446e-16 1.582468e-10
## 140 2.220446e-16 2.220446e-16 1.822660e-15
## 141 1.723626e-03 2.585469e-04 7.691971e-03
## 144 6.612362e-06 4.465474e-11 1.727118e-04
## 145 1.714537e-06 1.978584e-13 7.268882e-05
## 147 2.220446e-16 2.220446e-16 3.367864e-10
## 148 1.159310e-06 4.224519e-14 5.522442e-05
## 150 2.025694e-04 1.365469e-06 1.930345e-03
## 152 1.386316e-11 2.220446e-16 2.572739e-08
## 153 1.827117e-10 2.220446e-16 1.420127e-07
## 154 2.483492e-12 2.220446e-16 8.108317e-09
## 155 2.220446e-16 2.220446e-16 2.139821e-10
## 156 6.471417e-05 6.455280e-08 8.539619e-04
## 158 2.220446e-16 2.220446e-16 1.771537e-15
## 159 1.715462e-13 2.220446e-16 1.416523e-09
## 160 2.220446e-16 2.220446e-16 3.760026e-12
## 163 1.209410e-08 2.220446e-16 2.565775e-06
## 165 2.220446e-16 2.220446e-16 5.132848e-16
## 166 1.193759e-07 2.220446e-16 1.322111e-05
## 167 1.923017e-05 1.432272e-09 3.740335e-04
## 168 1.955610e-09 2.220446e-16 7.482823e-07
## 170 8.789743e-09 2.220446e-16 1.997092e-06
## 171 1.056595e-03 8.910836e-05 5.591604e-03
## 172 1.296358e-06 6.140682e-14 6.048446e-05
## 173 1.773635e-03 2.816789e-04 7.765272e-03
## 174 1.466513e-03 1.644356e-04 7.189116e-03
## 177 7.797777e-05 1.008999e-07 9.939507e-04
## 179 2.220446e-16 2.220446e-16 1.684655e-12
## 180 1.579200e-04 6.891840e-07 1.635691e-03
## 183 3.743746e-13 2.220446e-16 2.289031e-09
## 185 5.464741e-04 1.772126e-05 3.712285e-03
## 186 8.362679e-12 2.220446e-16 1.955656e-08
## 187 4.559914e-13 2.220446e-16 2.626427e-09
## 188 2.220446e-16 2.220446e-16 6.999974e-11
## 189 2.220446e-16 2.220446e-16 2.070321e-10
## 190 2.326855e-05 2.383578e-09 4.369662e-04
## 191 5.141178e-10 2.220446e-16 3.064074e-07
## 192 2.965905e-05 5.203554e-09 5.158448e-04
## 193 2.220446e-16 2.220446e-16 1.132190e-11
## 194 5.622548e-04 2.075706e-05 3.680994e-03
## 197 2.220446e-16 2.220446e-16 2.220446e-16
## 198 2.220446e-16 2.220446e-16 2.950286e-13
## 199 1.816024e-13 2.220446e-16 1.451592e-09
## 200 5.110459e-03 1.725580e-03 1.706068e-02
## 201 2.220446e-16 2.220446e-16 2.397014e-11
## 203 2.220446e-16 2.220446e-16 4.959688e-14
## 205 2.220446e-16 2.220446e-16 6.589720e-14
## 206 2.220446e-16 2.220446e-16 9.214788e-11
## 207 1.880366e-05 9.508327e-10 3.960759e-04
## 208 2.127206e-07 2.220446e-16 1.678987e-05
## 209 3.819611e-06 3.652248e-12 1.303836e-04
## 211 2.220446e-16 2.220446e-16 1.413839e-12
## 212 2.220446e-16 2.220446e-16 9.759619e-12
## 213 2.220446e-16 2.220446e-16 1.907851e-11
## 214 3.920842e-13 2.220446e-16 2.239946e-09
## 215 7.808573e-04 5.172864e-05 4.375164e-03
## 217 5.014992e-13 2.220446e-16 2.610644e-09
## 218 6.003980e-05 4.696893e-08 8.302907e-04
## 219 2.543099e-13 2.220446e-16 1.808420e-09
## 220 8.750816e-10 2.220446e-16 4.129484e-07
## 221 6.397080e-13 2.220446e-16 3.265214e-09
## 222 2.220446e-16 2.220446e-16 1.038279e-14
## 223 2.220446e-16 2.220446e-16 6.060215e-10
## 224 2.220446e-16 2.220446e-16 1.161734e-15
## 226 3.099889e-10 2.220446e-16 2.094095e-07
## 227 4.090708e-10 2.220446e-16 2.577290e-07
## 228 1.819931e-03 2.331832e-04 8.608295e-03
## 229 5.661764e-09 2.220446e-16 1.622207e-06
## 230 2.220446e-16 2.220446e-16 1.649392e-11
## 231 3.186384e-05 5.663474e-09 5.595689e-04
## 232 1.320953e-05 2.320171e-10 3.254548e-04
## 233 2.220446e-16 2.220446e-16 2.576676e-12
## 234 2.220446e-16 2.220446e-16 5.976223e-15
## 235 8.671616e-07 1.163557e-14 4.594824e-05
## 237 4.595098e-03 1.542822e-03 1.539350e-02
## 238 2.220446e-16 2.220446e-16 2.220446e-16
## 239 2.220446e-16 2.220446e-16 5.080541e-10
## 240 1.547359e-11 2.220446e-16 2.873067e-08
## 241 2.417770e-10 2.220446e-16 1.898716e-07
## 242 2.220446e-16 2.220446e-16 3.232180e-10
## 245 2.220446e-16 2.220446e-16 5.422582e-12
## 246 1.499827e-09 2.220446e-16 6.091559e-07
## 250 5.967602e-03 2.624721e-03 1.778279e-02
## 252 2.220446e-16 2.220446e-16 1.099050e-11
## 253 6.533937e-11 2.220446e-16 7.795234e-08
## 255 2.713310e-08 2.220446e-16 4.664209e-06
## 256 2.142374e-10 2.220446e-16 1.738299e-07
## 258 2.391171e-08 2.220446e-16 4.143149e-06
## 259 2.220446e-16 2.220446e-16 1.297868e-14
## 260 6.180454e-06 3.631885e-11 1.637316e-04
## 261 3.714659e-03 9.648521e-04 1.378575e-02
## 262 5.248088e-08 2.220446e-16 7.315492e-06
## 264 1.019415e-05 1.898734e-10 2.359649e-04
## 266 3.894488e-09 2.220446e-16 1.276747e-06
## 267 1.577613e-04 5.012537e-07 1.775853e-03
## 268 5.838716e-05 5.198613e-08 7.791498e-04
## 270 7.232816e-10 2.220446e-16 3.893433e-07
## 271 8.428557e-04 5.789857e-05 4.694854e-03
## 272 1.284891e-07 2.220446e-16 1.294719e-05
## 273 5.705521e-13 2.220446e-16 3.114709e-09
## 274 2.778960e-06 1.620330e-12 9.729982e-05
## 275 1.848776e-11 2.220446e-16 3.293709e-08
## 276 2.317551e-12 2.220446e-16 7.661620e-09
## 277 2.220446e-16 2.220446e-16 5.681533e-10
## 278 5.545457e-06 1.102828e-11 1.778182e-04
## 279 3.584743e-09 2.220446e-16 1.150390e-06
## 280 5.958616e-08 2.220446e-16 7.269954e-06
## 281 2.429141e-06 7.917735e-13 9.203206e-05
## 282 2.217223e-13 2.220446e-16 1.502458e-09
## 283 5.336262e-13 2.220446e-16 2.818209e-09
## 284 2.220446e-16 2.220446e-16 1.295883e-12
## 285 2.220446e-16 2.220446e-16 2.384345e-15
## 286 2.220446e-16 2.220446e-16 1.593853e-13
## 288 1.067606e-03 9.804411e-05 5.495545e-03
## 290 9.753463e-07 1.133968e-14 5.441787e-05
## 291 4.177819e-05 1.367668e-08 6.684951e-04
## 292 6.249918e-07 2.230145e-15 3.837998e-05
## 296 2.220446e-16 2.220446e-16 8.040320e-10
## 297 9.861878e-11 2.220446e-16 9.687439e-08
## 298 2.311734e-06 8.050785e-13 8.559584e-05
## 299 9.357287e-08 2.220446e-16 1.036730e-05
## 300 2.220446e-16 2.220446e-16 8.209454e-15
## 302 2.689135e-04 3.452664e-06 2.220662e-03
## 303 4.204298e-08 2.220446e-16 6.496955e-06
## 304 2.220446e-16 2.220446e-16 1.594223e-11
## 305 2.220446e-16 2.220446e-16 2.725139e-11
## 306 2.220446e-16 2.220446e-16 2.518565e-12
## 307 9.105819e-05 1.574744e-07 1.105736e-03
## 308 1.371200e-07 2.220446e-16 1.343215e-05
## 312 4.403931e-06 6.212241e-12 1.438214e-04
## 313 3.388123e-13 2.220446e-16 2.265094e-09
## 314 2.220446e-16 2.220446e-16 8.495615e-15
## 315 3.340755e-11 2.220446e-16 4.592899e-08
## 316 2.220446e-16 2.220446e-16 4.497967e-16
## 317 1.213042e-13 2.220446e-16 1.073358e-09
## 318 6.614624e-13 2.220446e-16 3.516415e-09
## 319 2.556611e-09 2.220446e-16 8.397166e-07
## 320 2.220446e-16 2.220446e-16 1.724489e-11
## 322 2.220446e-16 2.220446e-16 4.312607e-10
## 323 2.417789e-12 2.220446e-16 8.286998e-09
## 324 2.220446e-16 2.220446e-16 1.171444e-10
## 325 2.220446e-16 2.220446e-16 8.430751e-10
## 326 2.220446e-16 2.220446e-16 5.266767e-15
## 327 5.373551e-11 2.220446e-16 6.840720e-08
## 330 1.786171e-08 2.220446e-16 3.613377e-06
## 333 8.217888e-03 3.827662e-03 2.361093e-02
## 334 1.126103e-10 2.220446e-16 1.117009e-07
## 335 3.881832e-10 2.220446e-16 2.593163e-07
## 336 1.250602e-03 1.112307e-04 6.572852e-03
## 337 2.220446e-16 2.220446e-16 2.018149e-10
## 340 6.452868e-06 2.630544e-11 1.854391e-04
## 342 2.220446e-16 2.220446e-16 2.220446e-16
## 343 2.461014e-06 8.620649e-13 9.226982e-05
## 344 2.220446e-16 2.220446e-16 3.693025e-10
## 345 7.671291e-03 3.817543e-03 2.147187e-02
## 347 2.220446e-16 2.220446e-16 2.220446e-16
## 348 2.220446e-16 2.220446e-16 2.965879e-13
## 349 3.633427e-09 2.220446e-16 1.101754e-06
## 350 3.069175e-09 2.220446e-16 9.785560e-07
## 353 2.734520e-03 5.906927e-04 1.084386e-02
## 355 4.754388e-06 1.062030e-11 1.443709e-04
## 356 2.220446e-16 2.220446e-16 2.220446e-16
## 357 2.526235e-03 5.153660e-04 1.021564e-02
## 358 1.415227e-07 2.220446e-16 1.341106e-05
## 359 1.682447e-10 2.220446e-16 1.389984e-07
## 362 1.368931e-13 2.220446e-16 1.153044e-09
## 363 2.220446e-16 2.220446e-16 1.612252e-11
## 365 7.511464e-09 2.220446e-16 1.776693e-06
## 366 9.958438e-12 2.220446e-16 2.165269e-08
## 367 2.220446e-16 2.220446e-16 8.799995e-11
## 368 3.579817e-09 2.220446e-16 1.158385e-06
## 369 3.660915e-13 2.220446e-16 2.295973e-09
## 370 3.641880e-08 2.220446e-16 5.269800e-06
## 371 2.220446e-16 2.220446e-16 6.839896e-11
## 372 2.714053e-12 2.220446e-16 8.307925e-09
## 373 2.220446e-16 2.220446e-16 9.946444e-11
## 374 3.685177e-08 2.220446e-16 5.558077e-06
## 375 2.220446e-16 2.220446e-16 2.220446e-16
## 376 1.455106e-09 2.220446e-16 5.953202e-07
## 377 5.634074e-07 1.982427e-15 3.384391e-05
## 378 2.220446e-16 2.220446e-16 7.115624e-13
## 379 2.211925e-13 2.220446e-16 1.606314e-09
## 380 1.450911e-08 2.220446e-16 2.805671e-06
## 381 2.220446e-16 2.220446e-16 5.593975e-14
## 382 1.451917e-03 1.653893e-04 7.075126e-03
## 383 2.220446e-16 2.220446e-16 7.017955e-15
## 384 2.220446e-16 2.220446e-16 3.358408e-13
## 386 1.693806e-11 2.220446e-16 3.158428e-08
## 387 1.694296e-09 2.220446e-16 6.499061e-07
## 388 1.517424e-03 1.728260e-04 7.412090e-03
## 391 6.476935e-09 2.220446e-16 1.206735e-06
## 392 1.900264e-03 2.844758e-04 5.987220e-03
## 393 2.484481e-03 4.166666e-04 7.588104e-03
## 394 5.515276e-03 2.024508e-03 1.250069e-02
## 395 2.220446e-16 2.220446e-16 2.080755e-13
## 396 4.744635e-06 6.282530e-12 1.117934e-04
## 397 9.607628e-05 1.546215e-07 8.404658e-04
## 398 2.220446e-16 2.220446e-16 9.486555e-12
## 401 1.569587e-12 2.220446e-16 4.248841e-09
## 402 2.770158e-07 2.220446e-16 1.551697e-05
## 403 2.220446e-16 2.220446e-16 2.437787e-14
## 404 4.464949e-10 2.220446e-16 2.006302e-07
## 405 1.404285e-06 5.190145e-14 4.890172e-05
## 406 7.140592e-04 3.258974e-05 3.139779e-03
## 407 2.220446e-16 2.220446e-16 1.708022e-11
## 408 6.298243e-07 3.597050e-15 2.513194e-05
## 409 2.220446e-16 2.220446e-16 9.878436e-13
## 411 1.921064e-08 2.220446e-16 2.499651e-06
## 413 8.623385e-06 9.311021e-11 1.515996e-04
## 414 6.090910e-07 2.251086e-15 2.593652e-05
## 415 1.127220e-03 9.227100e-05 4.249832e-03
## 416 4.460663e-09 2.220446e-16 9.165930e-07
## 417 5.661024e-05 3.078049e-08 5.936669e-04
## 418 1.768682e-06 2.437662e-13 5.131597e-05
## 419 7.249040e-03 2.871177e-03 1.581128e-02
## 420 2.220446e-16 2.220446e-16 6.807556e-13
## 421 2.477775e-05 2.898828e-09 3.206923e-04
## 422 2.220446e-16 2.220446e-16 7.434757e-11
## 423 2.220446e-16 2.220446e-16 7.706672e-16
## 425 2.255472e-07 2.220446e-16 1.366829e-05
## 427 9.029292e-11 2.220446e-16 6.673691e-08
## 428 4.194978e-06 6.862196e-12 9.249336e-05
## 429 9.454154e-04 5.424111e-05 3.954527e-03
## 430 2.893779e-12 2.220446e-16 6.483701e-09
## 431 8.374919e-07 8.955553e-15 3.211179e-05
## 432 2.257987e-09 2.220446e-16 5.852839e-07
## 434 3.893540e-09 2.220446e-16 8.914014e-07
## 435 6.331960e-03 2.936350e-03 1.295607e-02
## 436 2.220446e-16 2.220446e-16 2.857636e-14
## 437 8.414918e-07 1.140941e-14 3.101754e-05
## 438 3.712670e-08 2.220446e-16 3.896574e-06
## 439 1.642983e-11 2.220446e-16 2.072271e-08
## 440 9.165706e-10 2.220446e-16 3.163280e-07
## 442 2.220446e-16 2.220446e-16 4.556123e-10
## 443 2.084265e-06 3.371616e-13 6.091237e-05
## 444 3.914931e-11 2.220446e-16 3.790496e-08
## 445 3.262119e-07 2.220446e-16 1.758642e-05
## 447 2.220446e-16 2.220446e-16 5.168843e-10
## 448 1.229674e-07 2.220446e-16 9.201337e-06
## 449 3.119872e-12 2.220446e-16 6.660283e-09
## 450 3.405394e-03 8.112736e-04 9.181281e-03
## 451 1.331088e-04 4.168936e-07 1.030711e-03
## 452 2.220446e-16 2.220446e-16 5.444897e-12
## 453 1.485913e-03 1.996959e-04 4.808383e-03
## 454 2.283535e-13 2.220446e-16 1.161205e-09
## 455 2.220446e-16 2.220446e-16 1.603401e-12
## 458 7.723219e-11 2.220446e-16 6.148578e-08
## 461 1.725208e-12 2.220446e-16 4.628578e-09
## 462 2.220446e-16 2.220446e-16 2.799140e-11
## 463 1.879043e-12 2.220446e-16 5.071683e-09
## 464 4.778238e-03 1.385085e-03 1.194635e-02
## 465 3.634936e-06 4.103501e-12 8.348669e-05
## 466 1.561680e-09 2.220446e-16 4.694514e-07
## 467 2.220446e-16 2.220446e-16 2.888034e-16
## 468 2.220446e-16 2.220446e-16 1.474772e-14
## 469 1.686198e-04 1.009114e-06 1.137015e-03
## 470 2.220446e-16 2.220446e-16 2.090534e-12
## 471 2.220446e-16 2.220446e-16 9.008172e-11
## 473 2.220446e-16 2.220446e-16 2.220446e-16
## 474 1.162572e-06 2.350281e-14 4.308449e-05
## 475 4.511808e-04 1.030931e-05 2.345243e-03
## 476 1.694282e-05 9.377628e-10 2.409039e-04
## 477 1.158959e-07 2.220446e-16 8.344088e-06
## 479 3.851081e-08 2.220446e-16 4.164941e-06
## 480 2.220446e-16 2.220446e-16 1.230019e-10
## 481 9.887020e-05 2.008754e-07 8.191723e-04
## 482 2.269068e-10 2.220446e-16 1.197456e-07
## 483 4.784403e-10 2.220446e-16 2.049227e-07
## 484 2.220446e-16 2.220446e-16 9.785904e-15
## 485 6.862828e-09 2.220446e-16 1.298120e-06
## 486 4.510284e-11 2.220446e-16 4.223936e-08
## 487 8.263198e-11 2.220446e-16 6.825973e-08
## 488 2.220446e-16 2.220446e-16 8.819674e-14
## 489 2.220446e-16 2.220446e-16 2.237778e-12
## 490 1.125152e-11 2.220446e-16 1.573823e-08
## 492 2.220446e-16 2.220446e-16 2.770075e-12
## 493 2.220446e-16 2.220446e-16 1.263110e-12
## 495 2.220446e-16 2.220446e-16 8.122650e-13
## 496 1.250627e-06 4.669715e-14 4.234827e-05
## 497 2.935528e-04 3.719461e-06 1.725793e-03
## 498 1.953516e-11 2.220446e-16 2.322373e-08
## 499 8.857202e-10 2.220446e-16 3.121972e-07
## 500 2.220446e-16 2.220446e-16 1.205348e-10
## 501 3.215663e-09 2.220446e-16 7.581166e-07
## 503 2.220446e-16 2.220446e-16 3.342227e-10
## 504 4.140917e-08 2.220446e-16 4.198433e-06
## 505 2.220446e-16 2.220446e-16 9.055462e-15
## 506 1.318836e-12 2.220446e-16 3.882870e-09
## 508 1.598244e-11 2.220446e-16 2.043255e-08
## 509 2.220446e-16 2.220446e-16 5.379125e-10
## 510 1.146522e-12 2.220446e-16 3.432291e-09
## 512 2.220446e-16 2.220446e-16 2.729488e-12
## 513 3.246052e-13 2.220446e-16 1.581158e-09
## 514 2.220446e-16 2.220446e-16 4.428135e-10
## 515 2.220446e-16 2.220446e-16 8.201113e-15
## 516 3.134313e-09 2.220446e-16 7.869755e-07
## 517 1.677503e-06 8.992843e-14 5.694316e-05
## 518 3.945519e-07 3.843582e-16 1.887492e-05
## 520 2.220446e-16 2.220446e-16 1.289386e-12
## 521 6.284978e-07 2.466334e-15 2.669230e-05
## 522 2.220446e-16 2.220446e-16 2.345635e-15
## 523 2.401190e-13 2.220446e-16 1.139062e-09
## 524 2.220446e-16 2.220446e-16 1.385111e-11
## 525 2.220446e-16 2.220446e-16 6.390326e-12
## 527 2.220446e-16 2.220446e-16 5.708255e-11
## 528 4.861092e-11 2.220446e-16 4.696244e-08
## 529 9.373821e-10 2.220446e-16 3.301692e-07
## 530 9.134676e-13 2.220446e-16 2.834459e-09
## 532 2.220446e-16 2.220446e-16 5.700926e-13
## 533 1.092099e-07 2.220446e-16 8.102587e-06
## 535 3.215784e-07 2.220446e-16 1.748921e-05
## 536 9.294255e-04 5.088111e-05 3.941635e-03
## 537 2.284511e-07 2.220446e-16 1.293141e-05
## 538 6.264969e-10 2.220446e-16 2.654713e-07
## 539 3.710070e-08 2.220446e-16 3.731212e-06
## 540 1.063361e-09 2.220446e-16 3.580713e-07
## 541 1.843013e-03 2.898537e-04 5.697879e-03
## 542 2.175033e-12 2.220446e-16 5.158842e-09
## 543 9.891053e-06 1.286734e-10 1.719770e-04
## 544 2.778415e-09 2.220446e-16 6.805014e-07
## 546 2.337119e-03 3.865061e-04 7.160025e-03
## 547 2.220446e-16 2.220446e-16 7.658573e-11
## 548 4.565816e-08 2.220446e-16 4.278765e-06
## 549 1.024692e-12 2.220446e-16 3.480113e-09
## 550 5.656303e-05 3.982491e-08 5.576394e-04
## 551 3.956308e-06 5.654196e-12 8.842991e-05
## 553 1.247469e-08 2.220446e-16 1.926349e-06
## 554 2.630356e-12 2.220446e-16 5.869618e-09
## 555 2.220446e-16 2.220446e-16 1.479879e-12
## 556 2.220446e-16 2.220446e-16 2.276053e-11
## 557 6.275440e-11 2.220446e-16 4.872935e-08
## 558 6.419277e-08 2.220446e-16 5.428701e-06
## 559 7.101998e-06 3.737451e-11 1.388101e-04
## 560 2.220446e-16 2.220446e-16 6.848163e-13
## 561 1.883446e-07 2.220446e-16 1.174985e-05
## 563 1.034642e-09 2.220446e-16 3.344077e-07
## 564 2.220446e-16 2.220446e-16 8.839747e-13
## 565 2.539469e-09 2.220446e-16 5.915835e-07
## 566 1.292424e-11 2.220446e-16 1.784461e-08
## 567 1.899933e-05 1.155874e-09 2.703721e-04
## 568 1.360965e-11 2.220446e-16 1.885085e-08
## 569 1.346836e-03 1.390140e-04 4.748807e-03
## 571 1.277561e-10 2.220446e-16 8.993456e-08
## 572 6.545709e-07 3.137650e-15 2.714181e-05
## 574 1.719311e-06 1.298906e-13 5.521346e-05
## 575 2.220446e-16 2.220446e-16 4.229853e-14
## 576 2.220446e-16 2.220446e-16 1.344069e-13
## 577 2.220446e-16 2.220446e-16 3.924279e-12
## 578 4.011375e-09 2.220446e-16 8.617600e-07
## 580 2.220446e-16 2.220446e-16 2.315155e-15
## 581 5.723705e-10 2.220446e-16 2.119747e-07
## 582 1.520088e-06 1.191894e-13 4.717785e-05
## 583 2.220446e-16 2.220446e-16 8.916114e-12
## 584 3.001911e-10 2.220446e-16 1.519346e-07
## 585 2.157955e-06 3.529543e-13 6.342280e-05
## 586 2.220446e-16 2.220446e-16 2.679049e-10
## 587 6.294001e-06 3.092333e-11 1.216983e-04
## 588 1.621279e-04 6.228321e-07 1.223407e-03
## 589 2.220446e-16 2.220446e-16 8.435248e-11
## 590 4.869363e-12 2.220446e-16 9.447050e-09
## 591 3.887761e-08 2.220446e-16 4.219852e-06
## 592 2.914870e-06 1.435729e-12 7.467103e-05
## 593 2.220446e-16 2.220446e-16 1.231754e-10
## 594 2.435881e-09 2.220446e-16 6.789379e-07
## 595 1.052991e-10 2.220446e-16 7.697297e-08
## 596 9.627915e-12 2.220446e-16 1.409431e-08
## 597 1.580855e-13 2.220446e-16 9.580764e-10
## 599 2.220446e-16 2.220446e-16 3.667970e-10
## 600 2.220446e-16 2.220446e-16 2.220446e-16
## 601 3.572667e-04 7.266435e-06 1.874575e-03
## 602 5.496816e-04 2.215537e-05 2.454196e-03
## 603 7.952874e-06 4.012760e-11 1.603736e-04
## 605 1.918694e-06 2.035607e-13 5.941513e-05
## 606 2.420571e-13 2.220446e-16 1.183793e-09
## 607 7.242075e-12 2.220446e-16 1.191954e-08
## 608 2.220446e-16 2.220446e-16 2.871733e-10
## 609 7.942417e-10 2.220446e-16 3.147100e-07
## 610 2.739178e-10 2.220446e-16 1.450394e-07
## 611 2.576430e-09 2.220446e-16 7.039769e-07
## 612 3.923813e-13 2.220446e-16 1.610851e-09
## 613 2.220446e-16 2.220446e-16 1.103476e-12
## 614 1.768515e-12 2.220446e-16 4.749625e-09
## 615 1.659398e-09 2.220446e-16 4.781885e-07
## 616 2.220446e-16 2.220446e-16 1.785129e-11
## 617 8.832319e-12 2.220446e-16 1.300787e-08
## 618 2.220446e-16 2.220446e-16 2.182105e-10
## 619 2.220446e-16 2.220446e-16 7.369109e-16
## 620 2.220446e-16 2.220446e-16 3.176092e-14
## 621 9.084634e-07 1.082437e-14 3.484395e-05
## 622 7.983198e-13 2.220446e-16 2.810029e-09
## 624 5.701294e-13 2.220446e-16 2.019950e-09
## 625 2.220446e-16 2.220446e-16 2.220446e-16
## 626 1.077137e-12 2.220446e-16 3.201299e-09
## 627 1.096891e-13 2.220446e-16 7.184023e-10
## 628 7.326798e-08 2.220446e-16 5.880298e-06
## 629 2.220446e-16 2.220446e-16 1.473703e-10
## 630 2.220446e-16 2.220446e-16 2.105154e-10
## 631 2.220446e-16 2.220446e-16 2.220446e-16
## 632 1.847491e-04 9.003472e-07 1.332300e-03
## 633 7.370829e-12 2.220446e-16 1.267044e-08
## 634 2.220446e-16 2.220446e-16 5.167306e-10
## 635 4.133452e-07 4.103459e-16 1.993397e-05
## 637 1.465030e-09 2.220446e-16 4.385614e-07
## 638 4.563186e-11 2.220446e-16 4.175669e-08
## 639 2.220446e-16 2.220446e-16 2.408810e-12
## 640 2.828607e-13 2.220446e-16 1.325653e-09
## 641 5.745730e-10 2.220446e-16 2.173293e-07
## 642 2.220446e-16 2.220446e-16 6.538032e-11
## 643 1.080845e-11 2.220446e-16 1.605101e-08
## 644 1.140921e-07 2.220446e-16 8.310908e-06
## 646 1.086414e-10 2.220446e-16 7.579812e-08
## 647 5.163867e-04 1.635369e-05 2.467195e-03
## 648 1.772298e-08 2.220446e-16 2.475933e-06
## 651 1.015607e-10 2.220446e-16 6.764129e-08
## 652 2.220446e-16 2.220446e-16 3.075084e-12
## 653 3.216869e-04 5.766658e-06 1.730703e-03
## 654 1.684373e-09 2.220446e-16 4.830064e-07
## 655 1.402623e-11 2.220446e-16 1.813054e-08
## 656 2.220446e-16 2.220446e-16 9.957073e-11
## 657 1.374539e-09 2.220446e-16 3.880338e-07
## 658 2.220446e-16 2.220446e-16 1.816312e-11
## 659 2.220446e-16 2.220446e-16 1.664978e-11
## 660 3.388317e-11 2.220446e-16 3.419079e-08
## 661 1.773073e-06 1.941309e-13 5.364835e-05
## 662 2.220446e-16 2.220446e-16 3.641813e-12
## 663 2.220446e-16 2.220446e-16 8.779728e-14
## 664 7.504757e-03 3.714277e-03 1.483692e-02
## 665 9.740286e-11 2.220446e-16 6.925748e-08
## 668 6.574347e-07 3.165968e-15 2.726728e-05
## 669 4.884939e-13 2.220446e-16 2.000751e-09
## 670 2.381634e-08 2.220446e-16 3.018749e-06
## 672 2.220446e-16 2.220446e-16 1.386590e-10
## 673 2.220446e-16 2.220446e-16 5.360688e-14
## 675 8.522520e-08 2.220446e-16 7.368833e-06
## 676 5.156945e-10 2.220446e-16 2.128427e-07
## 677 2.220446e-16 2.220446e-16 5.961954e-10
## 678 4.373223e-13 2.220446e-16 1.934836e-09
## 679 5.308541e-04 1.540597e-05 2.611322e-03
## 680 2.553240e-07 2.220446e-16 1.388545e-05
## 681 2.377695e-07 2.220446e-16 1.402658e-05
## 682 7.264213e-13 2.220446e-16 2.523535e-09
## 683 6.804719e-04 3.346911e-05 2.910460e-03
## 684 2.096001e-08 2.220446e-16 2.717885e-06
## 685 2.220446e-16 2.220446e-16 2.220446e-16
## 686 1.639187e-13 2.220446e-16 9.546255e-10
## 687 1.463939e-07 2.220446e-16 9.842818e-06
## 688 6.831606e-08 2.220446e-16 5.835931e-06
## 689 5.201403e-05 3.380813e-08 5.156100e-04
## 690 2.220446e-16 2.220446e-16 2.935227e-11
## 691 3.230503e-10 2.220446e-16 1.603166e-07
## 692 7.017060e-10 2.220446e-16 2.698312e-07
## 694 2.220446e-16 2.220446e-16 4.189626e-10
## 695 6.684135e-03 2.909480e-03 1.403241e-02
## 696 1.358462e-04 5.743024e-07 9.754918e-04
## 697 2.220446e-16 2.220446e-16 1.442590e-10
## 698 2.919643e-08 2.220446e-16 3.503755e-06
## 699 2.220446e-16 2.220446e-16 1.204220e-14
## 700 2.220446e-16 2.220446e-16 1.856541e-14
## 701 1.348769e-07 2.220446e-16 9.070081e-06
## 702 7.234235e-13 2.220446e-16 2.518616e-09
## 703 2.415305e-12 2.220446e-16 6.245247e-09
## 705 3.381381e-13 2.220446e-16 1.621165e-09
## 706 2.220446e-16 2.220446e-16 1.993050e-12
## 707 1.843330e-04 1.337966e-06 1.193793e-03
## 709 2.061114e-04 1.477018e-06 1.358938e-03
## 710 2.241322e-08 2.220446e-16 2.986378e-06
## 711 2.220446e-16 2.220446e-16 8.149849e-14
## 713 2.220446e-16 2.220446e-16 1.718778e-12
## 714 1.361723e-12 2.220446e-16 4.138514e-09
## 715 2.220446e-16 2.220446e-16 4.701100e-12
## 716 8.182899e-05 1.161686e-07 7.207623e-04
## 717 3.638699e-13 2.220446e-16 1.544635e-09
## 718 2.220446e-16 2.220446e-16 3.001628e-10
## 720 1.067691e-09 2.220446e-16 3.305311e-07
## 721 2.220446e-16 2.220446e-16 4.934177e-11
## 722 5.954935e-13 2.220446e-16 2.314106e-09
## 723 7.550896e-08 2.220446e-16 6.231795e-06
## 724 2.220446e-16 2.220446e-16 1.039797e-12
## 725 9.121971e-11 2.220446e-16 6.849310e-08
## 726 2.220446e-16 2.220446e-16 5.740348e-12
## 727 2.220446e-16 2.220446e-16 2.135994e-12
## 728 6.469419e-10 2.220446e-16 2.445141e-07
## 729 2.220446e-16 2.220446e-16 1.200770e-10
## 730 2.220446e-16 2.220446e-16 7.338279e-13
## 731 3.319230e-11 2.220446e-16 3.298763e-08
## 732 7.475606e-03 3.822164e-03 1.457054e-02
## 733 4.117251e-05 1.605035e-08 4.428682e-04
## 735 4.016139e-09 2.220446e-16 8.845595e-07
## 736 3.682677e-12 2.220446e-16 7.653777e-09
## 737 2.220446e-16 2.220446e-16 2.013463e-12
## 738 2.220446e-16 2.220446e-16 2.220446e-16
## 739 2.220446e-16 2.220446e-16 2.662279e-15
## 740 4.396436e-06 5.671172e-12 1.024594e-04
## 741 1.382994e-08 2.220446e-16 1.944662e-06
## 742 2.220446e-16 2.220446e-16 6.685970e-10
## 743 3.988875e-12 2.220446e-16 7.746964e-09
## 745 2.220446e-16 2.220446e-16 1.895559e-11
## 746 1.481536e-06 8.045442e-14 4.879624e-05
## 747 2.591657e-06 1.328266e-12 6.425575e-05
## 749 3.583774e-07 2.442636e-16 1.775020e-05
## 750 7.757063e-13 2.220446e-16 2.791553e-09
## 751 1.262015e-05 2.693015e-10 2.076905e-04
## 752 2.220446e-16 2.220446e-16 3.681487e-16
## 754 2.220446e-16 2.220446e-16 3.080287e-14
## 755 1.200182e-05 2.560211e-10 1.956242e-04
## 756 2.220446e-16 2.220446e-16 2.505307e-15
## 757 1.344137e-04 4.144648e-07 1.046480e-03
## 758 1.221912e-12 2.220446e-16 3.700388e-09
## 759 5.290630e-06 1.176295e-11 1.153663e-04
## 760 1.754663e-05 9.772708e-10 2.507817e-04
## 762 2.220446e-16 2.220446e-16 4.936403e-10
## 763 8.138400e-05 1.031916e-07 7.365542e-04
## 765 3.616997e-11 2.220446e-16 3.847422e-08
## 766 2.220446e-16 2.220446e-16 2.220446e-16
## 768 5.903997e-09 2.220446e-16 1.134033e-06
## 769 5.439077e-10 2.220446e-16 2.131948e-07
## 770 2.220446e-16 2.220446e-16 2.220446e-16
## 772 5.831696e-06 1.887855e-11 1.205542e-04
## 773 2.220446e-16 2.220446e-16 2.917070e-10
## 775 2.572772e-08 2.220446e-16 3.052340e-06
## 776 2.220446e-16 2.220446e-16 4.544482e-12
## 777 5.556322e-07 2.100274e-15 2.306835e-05
## 779 4.804859e-11 2.220446e-16 4.601570e-08
## 780 1.012667e-08 2.220446e-16 1.715835e-06
## 781 1.507452e-11 2.220446e-16 2.007973e-08
## 782 1.331697e-11 2.220446e-16 1.788738e-08
## 783 2.197662e-13 2.220446e-16 1.231836e-09Inference
- Initial values of logit and cloglog are similar and close, where as probit values are comparatively smaller that the other two.
- In the middle range, logit and probit take lead in turns but cloglog has higher probability.
- Lower tail: cloglog and logit have similar values, whereas probit has lower probability than the other two.
Step 15. Comparing the three models by showing their ROC curves in the same plot (see Step 12).
#probit
linpred.pr<-predict(lmod.s_ahsg_probit) # Linear predictor
predprob.pr<-predict(lmod.s_ahsg_probit,type ="response") #predicted probabilities
predout.pr<-(ifelse(predprob.pr<0.999,"no","yes")) ### Predicted outcomes using 0.5 as the threshold
FAAlog_R_wsa.pr<-data.frame(FAAlog,predprob.pr,predout.pr)
xtabs(~risky.landing+predout.pr,FAAlog_R_wsa.pr)## predout.pr
## risky.landing no yes
## 0 716 7
## 1 3 57threshF<-seq(0.01,0.999,0.01)
sensitivity.pr<-specificity.pr<-rep(NA,length(threshF))
for(j in seq (along=threshF)){
pp.pr<-ifelse(FAAlog_R_wsa.pr$predprob.pr< threshF[j],"no","yes")
xx.pr<-xtabs(~risky.landing+pp.pr,FAAlog_R_wsa.pr)
specificity.pr[j]<-xx.pr[1,1]/(xx.pr[1,1]+xx.pr[1,2])
sensitivity.pr[j]<-xx.pr[2,2]/(xx.pr[2,1]+xx.pr[2,2])
}
plot(1-specificity.pr,sensitivity.pr,type="l")
#cll
linpred.cll<-predict(lmod.s_ahsg_cloglog) # Linear predictor
predprob.cll<-predict(lmod.s_ahsg_cloglog,type ="response") #predicted probabilities
predout.cll<-(ifelse(predprob.cll<0.5,"no","yes")) ### Predicted outcomes using 0.5 as the threshold
FAAlog_R_wsa.cll<-data.frame(FAAlog,predprob.cll,predout.cll)
xtabs(~risky.landing+predout.cll,FAAlog_R_wsa.cll)## predout.cll
## risky.landing no yes
## 0 685 38
## 1 0 60threshF<-seq(0.01,0.999,0.01)
sensitivity.cll<-specificity.cll<-rep(NA,length(threshF))
for(j in seq (along=threshF)){
pp.cll<-ifelse(FAAlog_R_wsa.cll$predprob.cll< threshF[j],"no","yes")
xx.cll<-xtabs(~risky.landing+pp.cll,FAAlog_R_wsa.cll)
specificity.cll[j]<-xx.cll[1,1]/(xx.cll[1,1]+xx.cll[1,2])
sensitivity.cll[j]<-xx.cll[2,2]/(xx.cll[2,1]+xx.cll[2,2])
}
plot(1-specificity.cll,sensitivity.cll,col="Red",type="l")
plot(1-specificityF,sensitivityF,type="l")
lines(1-specificity.pr,sensitivity.pr,col="Blue")
lines(1-specificity.cll,sensitivity.cll,col="Red")
Step 16. Using each model to identify the top 5 risky landings. Do they point to the same flights?
dfc<-data.frame(predval)
dfc <- mutate(dfc, id = rownames(dfc))
FAA_U_RmvNA<- mutate(FAA_U_RmvNA, id = rownames(FAA_U_RmvNA))
head(arrange(dfc,logit),5)#5,7,8,10,11## logit probit cloglog id
## 1 2.220446e-16 2.220446e-16 8.322862e-11 5
## 2 2.220446e-16 2.220446e-16 8.654803e-13 7
## 3 2.220446e-16 2.220446e-16 2.117752e-12 8
## 4 2.220446e-16 2.220446e-16 5.111611e-10 10
## 5 2.220446e-16 2.220446e-16 1.792280e-11 11head(arrange(dfc,probit),5)#3,5,6,7,8## logit probit cloglog id
## 1 2.649443e-12 2.220446e-16 8.265371e-09 3
## 2 2.220446e-16 2.220446e-16 8.322862e-11 5
## 3 1.323259e-08 2.220446e-16 2.878625e-06 6
## 4 2.220446e-16 2.220446e-16 8.654803e-13 7
## 5 2.220446e-16 2.220446e-16 2.117752e-12 8head(arrange(dfc,cloglog),5)#89,105,197,238,342## logit probit cloglog id
## 1 2.220446e-16 2.220446e-16 2.220446e-16 89
## 2 2.220446e-16 2.220446e-16 2.220446e-16 105
## 3 2.220446e-16 2.220446e-16 2.220446e-16 197
## 4 2.220446e-16 2.220446e-16 2.220446e-16 238
## 5 2.220446e-16 2.220446e-16 2.220446e-16 342head(arrange(dfc,desc(logit)),5)#1## logit probit cloglog id
## 1 1 1 1 64
## 2 1 1 1 364
## 3 1 1 1 309
## 4 1 1 1 136
## 5 1 1 1 178head(arrange(dfc,desc(probit)),5)#1## logit probit cloglog id
## 1 1.0000000 1 1 19
## 2 0.9999999 1 1 29
## 3 1.0000000 1 1 30
## 4 1.0000000 1 1 56
## 5 1.0000000 1 1 64head(arrange(dfc,desc(cloglog)),5)#1## logit probit cloglog id
## 1 0.9996897 0.9999947 1 1
## 2 0.9998698 0.9999993 1 18
## 3 1.0000000 1.0000000 1 19
## 4 0.9999999 1.0000000 1 29
## 5 1.0000000 1.0000000 1 30head(arrange(dfc,desc(logit),desc(probit),desc(cloglog)),50)#1## logit probit cloglog id
## 1 1.0000000 1.0000000 1 64
## 2 1.0000000 1.0000000 1 364
## 3 1.0000000 1.0000000 1 309
## 4 1.0000000 1.0000000 1 136
## 5 1.0000000 1.0000000 1 178
## 6 1.0000000 1.0000000 1 410
## 7 1.0000000 1.0000000 1 389
## 8 1.0000000 1.0000000 1 181
## 9 1.0000000 1.0000000 1 56
## 10 1.0000000 1.0000000 1 385
## 11 1.0000000 1.0000000 1 645
## 12 1.0000000 1.0000000 1 19
## 13 1.0000000 1.0000000 1 269
## 14 1.0000000 1.0000000 1 341
## 15 1.0000000 1.0000000 1 244
## 16 1.0000000 1.0000000 1 30
## 17 1.0000000 1.0000000 1 265
## 18 1.0000000 1.0000000 1 329
## 19 1.0000000 1.0000000 1 753
## 20 1.0000000 1.0000000 1 161
## 21 0.9999999 1.0000000 1 764
## 22 0.9999999 1.0000000 1 29
## 23 0.9999999 1.0000000 1 771
## 24 0.9999997 1.0000000 1 204
## 25 0.9999997 1.0000000 1 671
## 26 0.9999997 1.0000000 1 251
## 27 0.9999997 1.0000000 1 91
## 28 0.9999997 1.0000000 1 623
## 29 0.9999996 1.0000000 1 162
## 30 0.9999993 1.0000000 1 310
## 31 0.9999984 1.0000000 1 142
## 32 0.9999969 1.0000000 1 331
## 33 0.9999951 1.0000000 1 552
## 34 0.9999926 1.0000000 1 767
## 35 0.9999917 1.0000000 1 97
## 36 0.9999871 1.0000000 1 311
## 37 0.9999827 1.0000000 1 511
## 38 0.9999806 1.0000000 1 748
## 39 0.9999590 1.0000000 1 124
## 40 0.9999496 1.0000000 1 115
## 41 0.9999416 0.9999999 1 570
## 42 0.9999392 0.9999999 1 169
## 43 0.9999265 0.9999999 1 666
## 44 0.9999098 0.9999998 1 649
## 45 0.9998698 0.9999993 1 18
## 46 0.9998310 0.9999988 1 263
## 47 0.9997773 0.9999973 1 562
## 48 0.9997540 0.9999966 1 339
## 49 0.9997015 0.9999939 1 545
## 50 0.9996951 0.9999942 1 243FAA_U_RmvNA%>%filter(id %in% c(64,364,309,136,178,410))#yes top 5 are same## # A tibble: 6 x 11
## aircraft duration no_pasg speed_ground speed_air height pitch distance
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 boeing 162. 72 129. 128. 33.9 4.14 5382.
## 2 boeing 209. 58 125. 126. 40.1 4.65 5147.
## 3 boeing 198. 68 127. 128. 23.8 2.99 5031.
## 4 boeing 155. 67 129. 128. 24.0 5.15 5058.
## 5 boeing 63.3 52 133. 133. 18.2 4.11 5343.
## 6 airbus 132. 60 131. 131. 28.3 3.66 4896.
## # ... with 3 more variables: long.landing <dbl>, risky.landing <dbl>,
## # id <chr>head(FAA_U_RmvNA,50)## # A tibble: 50 x 11
## aircraft duration no_pasg speed_ground speed_air height pitch distance
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 boeing 98.5 53 108. 109. 27.4 4.04 3370.
## 2 boeing 126. 69 102. 103. 27.8 4.12 2988.
## 3 boeing 112. 61 71.1 NA 18.6 4.43 1145.
## 4 boeing 197. 56 85.8 NA 30.7 3.88 1664.
## 5 boeing 90.1 70 59.9 NA 32.4 4.03 1050.
## 6 boeing 138. 55 75.0 NA 41.2 4.20 1627.
## 7 boeing 73.0 54 54.4 NA 24.0 3.84 805.
## 8 boeing 52.9 57 57.1 NA 19.4 4.64 574.
## 9 boeing 156. 61 85.4 NA 35.4 4.23 1699.
## 10 boeing 177. 56 61.8 NA 36.7 4.18 1138.
## # ... with 40 more rows, and 3 more variables: long.landing <dbl>,
## # risky.landing <dbl>, id <chr>If we observe the top 50 observations:
we notice that cloglog have 1 probability for more number of top line itmes Logit values tend to fall below 1 sooner that the values of probit.
This is very much inline to theoratical understanding of the concept of these functions.
5.2 Confidence Interval
Step 17. Use the probit model and hazard model to make prediction for the flight described in Step 13. Report the predicted probability as well as its 95% confidence interval. Compare the results with that from Step 13.
#probit model
predict(lmod.s_ahsg_probit,newdata=new.ind,type ="link") ### Linear predictor## 1
## 9.321099predict(lmod.s_ahsg_probit,newdata= new.ind,type ="response") ### Probability## 1
## 1predict(lmod.s_ahsg_probit,newdata=new.ind,type ="link",se=T) ### Linear predictor## $fit
## 1
## 9.321099
##
## $se.fit
## [1] 1.70394
##
## $residual.scale
## [1] 1round(ilogit(c(9.321099+1.96*1.70394,9.321099+1.96*1.70394)),5 ) ### Confidence interval## [1] 1 1predict(lmod.s_ahsg_probit,newdata= new.ind,type ="response",se=T) ### Probability## $fit
## 1
## 1
##
## $se.fit
## 1
## 3.783508e-16
##
## $residual.scale
## [1] 1round(c(1+1.96*3.783508e-16,1-1.96*3.783508e-16), 5 ) ### Confidence interval## [1] 1 1#cll model
predict(lmod.s_ahsg_cloglog,newdata=new.ind,type ="link") ### Linear predictor## 1
## 11.10962predict(lmod.s_ahsg_cloglog,newdata= new.ind,type ="response") ### Probability## 1
## 1### Predict the value with its standard error
predict(lmod.s_ahsg_cloglog,newdata=new.ind,type ="link",se=T) ### Linear predictor## $fit
## 1
## 11.10962
##
## $se.fit
## [1] 1.958174
##
## $residual.scale
## [1] 1round(ilogit(c( 11.10962+1.96*1.958174,11.10962+1.96*1.958174)),5 ) ### Confidence interval## [1] 1 1predict(lmod.s_ahsg_cloglog,newdata= new.ind,type ="response",se=T) ### Probability## $fit
## 1
## 1
##
## $se.fit
## 1
## 4.348019e-16
##
## $residual.scale
## [1] 1round(c(1+1.96*4.348019e-16,1-1.96*4.348019e-16), 5 ) ### Confidence interval## [1] 1 1We observe that confidence Interval is 1,1 which is an indicator of better prediction ability of our selected model
6.0 Conclusion
Is logit better than probit, or vice versa? Both methods will yield similar (though not identical) inferences. Logit - also known as logistic regression - is more popular in health sciences like epidemiology partly because coefficients can be interpreted in terms of odds ratios. Probit models can be generalized to account for non-constant error variances in more advanced econometric settings (known as heteroskedastic probit models) and hence are used in some contexts by economists and political scientists. If these more advanced applications are not of relevance, than it does not matter which method you choose to go with.