BMishra Dissertation BWS Data Analysis
Bijesh Mishra
9/23/2021
- This document is developed by Bijesh MISHRA for dissertation based on the tutorial from a paper Aizaki and Fogarty, 2019. An R package and tutorial for case 2 best–worst scaling, Journal of Choice Modelling. DOI: https://doi.org/10.1016/j.jocm.2019.100171. This document uses “support.BWS2” package in R.
support.bws2 package: Tool for Case 2 best-worst scaling:
- This pacakge has three basic functions for case 2 (profile case) best-worst scaling:
- Convert orthogonal main-effect design into questions
- Create a dataset suitable for analysis
- Canculate count-based scores
Construct a profile using an orthogonal main effect design (OMED). A profile has K attributes and each attribute has \(L_K\) levels. If all the attributes have the same number of levels L, \(L^k\) omed is used to construct profile. If a profile has four attributes (A, B, C, D) each with three levels (A1, A2, A3, B1, B2, B3, C1, C2, C3, D1, D2, D3), \(3^4\) OMED are generated.
Example has four attributes and three levels of each attribute in the format of “A1: Attribute A Level 1.”:
A) Deer Sanctuary: A1: No Deer Sanctuary No Food Plot, A2: Deer Sanctuary without Food Plot, A3: Deer Sanctuary with Food Plot
B) Tree Canopy Cover: B1: Closed Canopy, B2: Moderate Canopy, B3: Open Canopy
C) Number of Deer Observed: C1: 6, C2: 10, C3: 16
D) Lease Fee ($): D1: USD 6, D2: USD 10, D3: USD 16
Total 18 BWS cases were constructed and divided into three different survey sets. Each respondent only receives 6 out of 18 sets.
Please note that your best and worst variables from each choise set may need to be named as “B#” and “W#” where, # is a number. Eg: B1 and W1 are Best and Worst levels and attributes respectively selected from choice set 1. Otherwise this package will give you a problem.
rm(list = ls()) # Clean environment.
options(digits = 4) # Print up to 4 digits.
setwd("/Users/bmishra/Dropbox/OSU/PhD/Dissertation/Survey/Data/FinalData/Objective2/R")
# Load Libraries:
library(support.BWS2, warn.conflicts = F, quietly = T)
library(dplyr, warn.conflicts = F, quietly = T)
library(survival, warn.conflicts = F, quietly = T) # clogit.
library(DoE.base, warn.conflicts = F, quietly = T)## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.designlibrary(openxlsx, warn.conflicts = F, quietly = T) # Write in Excel
library(readxl, warn.conflicts = F, quietly = T)
library(stargazer, warn.conflicts = F, quietly = T) #Combined Model Outputs.##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerWork Flow:
Step 1: Define choice situation: See Choice sets below and Design Matrix.
# Step 2: Creation of Design Matrix, profiles and Choice sets:
attr.lev.s = list(
sanctuary = c("nsnfp", "snpf", "sfp"),
canopy = c("open", "moderate", "close"),
deer = c("deer1", "deer6", "deer10"),
lease = c("les6", "les10", "les16"))
# Define base levels of each attribute:
base.lev.s = list(sanctuary = c("nsnfp"),
canopy = c("open"),
deer = c("deer1"),
lease = c("les6"))
# Step 3: Load data from the package. This is data collected using survey.
# Load data:
obj2.data = read_excel("/Users/bmishra/Dropbox/OSU/PhD/Dissertation/Survey/Data/FinalData/SurveyData02192021.xlsx", sheet = "BWData")
# Complete Orthogonal Design Matrix for all 18 Sets:
des.s = cbind(
#c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18)
#c(B1, B2, B3, B4, B5, B6, B1, B2, B3, B4, B5, B6, B1, B2, B3, B4, B5, B6)
c(2, 3, 1, 3, 2, 1, 3, 3, 1, 2, 1, 2, 1, 2, 1, 3, 2, 3),
c(1, 2, 1, 3, 3, 2, 1, 2, 1, 3, 3, 2, 2, 2, 3, 1, 1, 3),
c(1, 1, 2, 2, 3, 3, 2, 3, 3, 2, 1, 1, 2, 2, 1, 1, 3, 3),
c(1, 3, 1, 2, 3, 2, 2, 1, 3, 3, 1, 2, 3, 1, 2, 3, 2, 1))
# BWS Survey Set 1:
# Orthogonal Design Matrix for Set 1:
des.s1 = cbind(
#c (1, 2, 3, 4, 5, 6)
#c(B1, B2, B3, B4, B5, B6)
c(2, 3, 1, 3, 2, 1),
c(1, 2, 1, 3, 3, 2),
c(1, 1, 2, 2, 3, 3),
c(1, 3, 1, 2, 3, 2))
obj2.set1 = subset(obj2.data, SurveySet == 1)
obj2.set1 = cbind(obs = seq(1:length(obj2.set1$SurveySet)),
obj2.set1)
# Define response variable:
response.vars.s1 = colnames(obj2.set1[16:27])
# Develop questionnaire sets (Set 1):
bws2.questionnaire(choice.sets = des.s1,
attribute.levels = attr.lev.s,
position = "center")##
## Q1
## Best Attribute Worst
## [ ] snpf [ ]
## [ ] open [ ]
## [ ] deer1 [ ]
## [ ] les6 [ ]
##
## Q2
## Best Attribute Worst
## [ ] sfp [ ]
## [ ] moderate [ ]
## [ ] deer1 [ ]
## [ ] les16 [ ]
##
## Q3
## Best Attribute Worst
## [ ] nsnfp [ ]
## [ ] open [ ]
## [ ] deer6 [ ]
## [ ] les6 [ ]
##
## Q4
## Best Attribute Worst
## [ ] sfp [ ]
## [ ] close [ ]
## [ ] deer6 [ ]
## [ ] les10 [ ]
##
## Q5
## Best Attribute Worst
## [ ] snpf [ ]
## [ ] close [ ]
## [ ] deer10 [ ]
## [ ] les16 [ ]
##
## Q6
## Best Attribute Worst
## [ ] nsnfp [ ]
## [ ] moderate [ ]
## [ ] deer10 [ ]
## [ ] les10 [ ]# BWS Survey Set 2:
# Orthogonal Design Matrix for Set 2:
des.s2 = cbind(
#c (7, 8, 9, 10, 11, 12)
#c(B1, B2, B3, B4, B5, B6)
c(3, 3, 1, 2, 1, 2),
c(1, 2, 1, 3, 3, 2),
c(2, 3, 3, 2, 1, 1),
c(2, 1, 3, 3, 1, 2))
obj2.set2 = subset(obj2.data, SurveySet == 1)
obj2.set2 = cbind(obs = seq(1:length(obj2.set2$SurveySet)),
obj2.set2)
# Define response variable:
response.vars.s2 = colnames(obj2.set2[16:27])
# Develop questionnaire sets (Set 2):
bws2.questionnaire(choice.sets = des.s2,
attribute.levels = attr.lev.s,
position = "center")##
## Q1
## Best Attribute Worst
## [ ] sfp [ ]
## [ ] open [ ]
## [ ] deer6 [ ]
## [ ] les10 [ ]
##
## Q2
## Best Attribute Worst
## [ ] sfp [ ]
## [ ] moderate [ ]
## [ ] deer10 [ ]
## [ ] les6 [ ]
##
## Q3
## Best Attribute Worst
## [ ] nsnfp [ ]
## [ ] open [ ]
## [ ] deer10 [ ]
## [ ] les16 [ ]
##
## Q4
## Best Attribute Worst
## [ ] snpf [ ]
## [ ] close [ ]
## [ ] deer6 [ ]
## [ ] les16 [ ]
##
## Q5
## Best Attribute Worst
## [ ] nsnfp [ ]
## [ ] close [ ]
## [ ] deer1 [ ]
## [ ] les6 [ ]
##
## Q6
## Best Attribute Worst
## [ ] snpf [ ]
## [ ] moderate [ ]
## [ ] deer1 [ ]
## [ ] les10 [ ]# BWS Survey Set 3:
# Orthogonal Design Matrix for Set 3:
des.s3 = cbind(
#c(13, 14, 15, 16, 17, 18)
#c(B1, B2, B3, B4, B5, B6)
c(1, 2, 1, 3, 2, 3),
c(2, 2, 3, 1, 1, 3),
c(2, 2, 1, 1, 3, 3),
c(3, 1, 2, 3, 2, 1))
obj2.set3 = subset(obj2.data, SurveySet == 1)
obj2.set3 = cbind(obs = seq(1:length(obj2.set3$SurveySet)),
obj2.set3)
# Define response variable:
response.vars.s3 = colnames(obj2.set3[16:27])
# Develop questionnaire sets (Set 3):
bws2.questionnaire(choice.sets = des.s3,
attribute.levels = attr.lev.s,
position = "center")##
## Q1
## Best Attribute Worst
## [ ] nsnfp [ ]
## [ ] moderate [ ]
## [ ] deer6 [ ]
## [ ] les16 [ ]
##
## Q2
## Best Attribute Worst
## [ ] snpf [ ]
## [ ] moderate [ ]
## [ ] deer6 [ ]
## [ ] les6 [ ]
##
## Q3
## Best Attribute Worst
## [ ] nsnfp [ ]
## [ ] close [ ]
## [ ] deer1 [ ]
## [ ] les10 [ ]
##
## Q4
## Best Attribute Worst
## [ ] sfp [ ]
## [ ] open [ ]
## [ ] deer1 [ ]
## [ ] les16 [ ]
##
## Q5
## Best Attribute Worst
## [ ] snpf [ ]
## [ ] open [ ]
## [ ] deer10 [ ]
## [ ] les10 [ ]
##
## Q6
## Best Attribute Worst
## [ ] sfp [ ]
## [ ] close [ ]
## [ ] deer10 [ ]
## [ ] les6 [ ]Step 4, 5 and 6: Dataset for Models:
Data formatting for the binary part of case 2 (Yes/No) is not available in this document now and is work on progress.
The datasets for the paired model and the marginal model were created using bws2.dataset() and then stored in the objects model1.data to model9.data.
Pay attention to formatting of STR variable in each dataset. This variable differs for paired, marginal and marginal sequencing model.
Note: Information below for paired, marginal, and marginal sequential models need to be verified. I have only done initial verification. So, I request to verify information using your own data before believing what I have written here. This information may change.
For paired model, all 12 pairs receives same value for STR variable. For Marginal model, four best choices receives one value and rest 4 worst value receives different value for each set. In marginal sequence models, first 4 attributes have same STR values and rest 3 worst attributes have same value but different than previous four attributes.
Dataset for Model 1
# Model 1: Response variable: 1 for selected pair of best and worse levels, 0 otherwise.
# Model 1 (Paired Model) Data:
# Survey Set 1:
model1.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
model = "paired")
# Survey Set 2:
model1.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
model = "paired")
# Survey Set 3:
model1.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
model = "paired")
# Model 1 Survey Set Combined:
model1.data = rbind(model1.set1, model1.set2, model1.set3)
# View(model1.data)Dataset for Model 2
# Model 2:
# Response variable (Same as Model 1): 1 for selected pair of best and worse levels, 0 otherwise.
# Model 2 (Paired Model) Data:
# Survey Set 1:
model2.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = TRUE,
model = "paired")
# Survey Set 2:
model2.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = TRUE,
model = "paired")
# Survey Set 3:
model2.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = TRUE,
model = "paired")
# Model 2 Survey Set Combined:
model2.data = rbind(model2.set1, model2.set2, model2.set3)
# View(model2.data)Dataset for Model 3
# Model 3:
# Response variable: 1 if a level is selected as best or worse, 0 otherwise.
# Model 3 (Marginal Model) Data:
# Survey Set 1:
model3.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
model = "marginal")
# Survey Set 2:
model3.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
model = "marginal")
# Survey Set 3:
model3.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
model = "marginal")
# Model 3 Survey Set Combined:
model3.data = rbind(model3.set1, model3.set2, model3.set3)
# View(model3.data)Dataset for Model 4, 5, 6
# Model 4, 5, 6:
# Response variable (Same as Model 3): 1 if a level is selected as best or worse, 0 otherwise.
# Model 4, 5, 6 (Marginal Model) Data:
# Survey Set 1:
model456.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = FALSE,
base.level = base.lev.s,
model = "marginal")
# Survey Set 2:
model456.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = FALSE,
base.level = base.lev.s,
model = "marginal")
# Survey Set 3:
model456.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = FALSE,
base.level = base.lev.s,
model = "marginal")
# Model 4, 5, 6 Survey Set Combined:
model456.data = rbind(model456.set1, model456.set2, model456.set3)
# View(model456.data)Dataset for Model 7
# Model 7:
# Response variable (Same as Model 3): 1 if a level is selected as best or worse, 0 otherwise.
# Model 7 (Marginal Sequential) Data:
# Survey Set 1:
model7.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = FALSE,
model = "sequential")
# Survey Set 2:
model7.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = FALSE,
model = "sequential")
# Survey Set 3:
model7.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = FALSE,
model = "sequential")
# Model 7 Survey Set Combined:
model7.data = rbind(model7.set1, model7.set2, model7.set3)
# View(model7.data)Dataset for Model 8
# Model 8:
# Response variable (Same as Model 3): 1 if a level is selected as best or worse, 0 otherwise.
# Model 8 (Marginal Sequential) Data:
# Survey Set 1:
model8.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
base.attribute = "canopy",
model = "sequential")
# Survey Set 2:
model8.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
base.attribute = "canopy",
model = "sequential")
# Survey Set 3:
model8.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = TRUE,
base.level = base.lev.s,
base.attribute = "canopy",
model = "sequential")
# Model 8 Survey Set Combined:
model8.data = rbind(model8.set1, model8.set2, model8.set3)
# View(model8.data)Dataset for Model 9
# Model 9:
# Response variable (Same as Model 3): 1 if a level is selected as best or worse, 0 otherwise.
# Model 9 (Marginal Sequential) Data:
# Survey Set 1:
model9.set1 = bws2.dataset(
data = obj2.set1,
id = "obs",
response = response.vars.s1,
choice.sets = des.s1,
attribute.levels = attr.lev.s,
reverse = FALSE,
base.level = base.lev.s,
model = "sequential")
# Survey Set 2:
model9.set2 = bws2.dataset(
data = obj2.set2,
id = "obs",
response = response.vars.s2,
choice.sets = des.s2,
attribute.levels = attr.lev.s,
reverse = FALSE,
base.level = base.lev.s,
model = "sequential")
# Survey Set 3:
model9.set3 = bws2.dataset(
data = obj2.set3,
id = "obs",
response = response.vars.s3,
choice.sets = des.s3,
attribute.levels = attr.lev.s,
reverse = FALSE,
base.level = base.lev.s,
model = "sequential")
# Model 9 Survey Set Combined:
model9.data = rbind(model9.set1, model9.set2, model9.set3)
# View(model9.data)Output data files to excel:
# Understanding Data for Various Models:
# View DataSets for Various Models:
# View(model1.data)
# View(model2.data)
# View(model3.data)
# View(model456.data)
# View(model7.data)
# View(model8.data)
# View(model9.data)
# Design Matrix of Datasets:
# model1.data[1:12,] # Paired
# model2.data[1:12,] # Paired
# model3.data[1:8,] # Marginal
# model456.data[1:8,] #Marginal
# model7.data[1:7,] #Marginal Sequential
# model8.data[1:7,] # Marginal Sequential
# model9.data[1:7,] # Marginal Sequential
# # Export Data in Permanent Storage Format (Excel and CSV):
# write.xlsx(model1.data, file = "Model1Data.xlsx", overwrite = T)
# write.xlsx(model2.data, file = "Model2Data.xlsx", overwrite = T)
# write.xlsx(model3.data, file = "Model3Data.xlsx", overwrite = T)
# write.xlsx(model456.data, file = "Model456Data.xlsx", overwrite = T)
# write.xlsx(model7.data, file = "Model7Data.xlsx", overwrite = T)
# write.xlsx(model8.data, file = "Model8Data.xlsx", overwrite = T)
# write.xlsx(model9.data, file = "Model9Data.xlsx", overwrite = T)
#
# write.csv(model1.data, file = "Model1Data.csv")
# write.csv(model2.data, file = "Model2Data.csv")
# write.csv(model3.data, file = "Model3Data.csv")
# write.csv(model456.data, file = "Model456Data.csv")
# write.csv(model7.data, file = "Model7Data.csv")
# write.csv(model8.data, file = "Model8Data.csv")
# write.csv(model9.data, file = "Model9Data.csv")Step 7: Analysis
Summary of Model 1:
# BWS Scores:
scores1 = bws2.count(model1.data)
dim(scores1)## [1] 264 155sum(scores1, "level")barplot(scores1, "bw", "level")
# Agregate Best and Worst Scores:
apply(model1.data[model1.data$RES == 1, c("BEST.LV", "WORST.LV")], 2, table)## BEST.LV WORST.LV
## close 32 72
## deer1 230 82
## deer10 307 40
## deer6 153 124
## les10 29 194
## les16 31 234
## les6 30 223
## moderate 103 22
## nsnfp 76 242
## open 66 47
## sfp 213 70
## snpf 164 84# # # Best Worst Score According to Group of Respondents:
# sum(scores1[obj2.data$F4Gender == 0, ], "level") # Male BW Score
# sum(scores1[obj2.data$F4Gender == 1, ], "level") #Female BW Score
#
# # BarPlot for Respondents in their 20s and those in their 50s:
# barplot(scores1[obj2.data$F4Gender == 0, ], "bw", "level") # Male BW Score
# barplot(scores1[obj2.data$F4Gender == 1, ], "bw", "level") # Female BW Score
# Conditional Logit Model:
model1.mf = clogit(RES ~ sanctuary + deer + lease +
snpf + sfp + moderate + close +
deer6 + deer10 + les10 + les16 + strata(STR),
data = model1.data)
model1.mf## Call:
## clogit(RES ~ sanctuary + deer + lease + snpf + sfp + moderate +
## close + deer6 + deer10 + les10 + les16 + strata(STR), data = model1.data)
##
## coef exp(coef) se(coef) z p
## sanctuary 1e-04 1e+00 5e-02 0 1.0
## deer 5e-01 2e+00 5e-02 10 <2e-16
## lease -8e-01 4e-01 5e-02 -15 <2e-16
## snpf 3e-01 1e+00 6e-02 6 2e-08
## sfp 6e-01 2e+00 6e-02 11 <2e-16
## moderate 3e-01 1e+00 6e-02 5 9e-08
## close -3e-01 7e-01 6e-02 -5 9e-08
## deer6 -6e-01 6e-01 6e-02 -10 <2e-16
## deer10 6e-01 2e+00 6e-02 10 <2e-16
## les10 8e-02 1e+00 6e-02 1 0.2
## les16 -5e-02 9e-01 5e-02 -1 0.3
##
## Likelihood ratio test=1042 on 11 df, p=<2e-16
## n= 17712, number of events= 1434
## (19656 observations deleted due to missingness)summary(model1.mf)## Call:
## coxph(formula = Surv(rep(1, 37368L), RES) ~ sanctuary + deer +
## lease + snpf + sfp + moderate + close + deer6 + deer10 +
## les10 + les16 + strata(STR), data = model1.data, method = "exact")
##
## n= 17712, number of events= 1434
## (19656 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sanctuary 0.000137 1.000137 0.051404 0.00 1.00
## deer 0.489609 1.631678 0.051746 9.46 < 2e-16 ***
## lease -0.803850 0.447602 0.053011 -15.16 < 2e-16 ***
## snpf 0.332416 1.394332 0.058885 5.65 1.7e-08 ***
## sfp 0.633261 1.883744 0.058405 10.84 < 2e-16 ***
## moderate 0.317280 1.373387 0.059390 5.34 9.2e-08 ***
## close -0.319236 0.726704 0.059789 -5.34 9.3e-08 ***
## deer6 -0.596157 0.550925 0.059722 -9.98 < 2e-16 ***
## deer10 0.568973 1.766452 0.055802 10.20 < 2e-16 ***
## les10 0.075544 1.078470 0.056281 1.34 0.18
## les16 -0.052656 0.948706 0.054895 -0.96 0.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sanctuary 1.000 1.000 0.904 1.106
## deer 1.632 0.613 1.474 1.806
## lease 0.448 2.234 0.403 0.497
## snpf 1.394 0.717 1.242 1.565
## sfp 1.884 0.531 1.680 2.112
## moderate 1.373 0.728 1.222 1.543
## close 0.727 1.376 0.646 0.817
## deer6 0.551 1.815 0.490 0.619
## deer10 1.766 0.566 1.583 1.971
## les10 1.078 0.927 0.966 1.204
## les16 0.949 1.054 0.852 1.056
##
## Concordance= 0.742 (se = 0.007 )
## Likelihood ratio test= 1042 on 11 df, p=<2e-16
## Wald test = 838 on 11 df, p=<2e-16
## Score (logrank) test = 971 on 11 df, p=<2e-16# The attribute-level variables are effect-coded ones, and thus the coefficient of the base level in each attribute can be calculated using:
# Define base levels of each attribute:
b.m1 = coef(model1.mf)
(nsnfp = -sum(b.m1[4:5]))## [1] -0.9657names(nsnfp) = "NSNFP"
(open = -sum(b.m1[6:7]))## [1] 0.001956names(open) = "Open"
(deer1 = -sum(b.m1[8:9]))## [1] 0.02718names(deer1) = "1 Deer"
(les6 = -sum(b.m1[10:11]))## [1] -0.02289names(les6) = "$6 Lease"
canopy = 0
names(canopy) = "Canopy"
model1 = c(b.m1[1:2], canopy, b.m1[3], b.m1[4:5],
nsnfp, b.m1[6:7], open, b.m1[8:9],
deer1, b.m1[10:11], les6)
model1## sanctuary deer Canopy lease snpf sfp NSNFP moderate
## 0.000137 0.489609 0.000000 -0.803850 0.332416 0.633261 -0.965677 0.317280
## close Open deer6 deer10 1 Deer les10 les16 $6 Lease
## -0.319236 0.001956 -0.596157 0.568973 0.027185 0.075544 -0.052656 -0.022888barplot(model1)
Summary of Model 2:
Summary of Model 3:
Summary of Model 4:
Summary of Model 5:
Summary of Model 6:
Summary of Model 7:
Summary of Model 8:
Summary of Model 9:
# As mentioned in Flynn et al. (2008), the results from the paired model
# are similar to those from the marginal model: the correlation coefficient
# for the two results is calculated as follows:
# pairs(cbind(model1, model2, model3,
# model4, model5, model6,
# model7, model8, model9))
# cor(cbind(model1, model2, model3,
# model4, model5, model6,
# model7, model8, model9))
# # lower.tri(a, diag = T)
# boxplot(cbind(model1, model2, model3,
# model4, model5, model6,
# model7, model8, model9))
# barplot(cbind(model1, model2, model3,
# model4, model5, model6,
# model7, model8, model9))