require(lavaan)
 bifactor <- "
 general.factor =~ Easy_Reservation + Preferred_Seats
  + Flight_Options + Ticket_Prices
 + Seat_Comfort + Seat_Roominess
 + Overhead_Storage + Clean_Aircraft
 + Courtesy + Friendliness + Helpfulness + Service
 
 ticketing =~ Easy_Reservation + Preferred_Seats + Flight_Options + Ticket_Prices
 aircraft  =~ Seat_Comfort + Seat_Roominess + Overhead_Storage + Clean_Aircraft
 service   =~ Courtesy + Friendliness + Helpfulness + Service"
 
 modelfit <- cfa(bifactor, data=ratings[,1:12], orthogonal=T)
 summary(modelfit, fit.measures = TRUE, standardize = T)

 bifactor <- "
 general =~ Easy_Reservation + Preferred_Seats 
  + Flight_Options + Ticket_Prices
 + Seat_Comfort + Seat_Roominess + Overhead_Storage 
 + Clean_Aircraft + Courtesy + Friendliness
 + Helpfulness + Service
 ticketing =~ Easy_Reservation + Preferred_Seats + Flight_Options + Ticket_Prices
 aircraft  =~ Seat_Comfort + Seat_Roominess + Overhead_Storage + Clean_Aircraft
 service   =~ Courtesy + Friendliness + Helpfulness + Service
 Satisfaction ~ general + ticketing + aircraft + service
 "
 satfit <- sem(bifactor, data=ratings[,c(1:12,13)], orthogonal=T)
 summary(satfit, fit.measures = T, standardize = T)
 inspect(satfit, "rsquare")

 library(OpenMx)
 # make up some helpful lists of items we will need
 serviceVars     = c("Courtesy", "Friendliness", "Helpfulness", "Service")
 aircraftVars    = c("Seat_Comfort", "Seat_Roominess", "Overhead_Storage", "Clean_Aircraft")
 ticketingVars   = c("Easy_Reservation", "Preferred_Seats", "Flight_Options", "Ticket_Prices")
 allMeasuredVars = c(serviceVars, aircraftVars, ticketingVars)
 allLatentVars   = c("general_factor", "ticketing", "aircraft","service")

 bifactor <- mxModel("bifactor", type="RAM",
 manifestVars = allMeasuredVars,
 latentVars   = allLatentVars,
 mxPath(from = allLatentVars  , arrows = 2, free = F, values = 1), # latents fixed@1
 mxPath(from = allMeasuredVars, arrows = 2), # manifest residuals
 # Factor loadings
 mxPath(from = "general_factor", to = allMeasuredVars),
 mxPath(from = "ticketing"     , to = ticketingVars),
 mxPath(from = "aircraft"      , to = aircraftVars),
 mxPath(from = "service"       , to = serviceVars),
 mxData(cov(ratings[,allMeasuredVars]), type = "cov", numObs = nrow(ratings))
 )
 bifactor= mxRun(bifactor)
 tmp = summary(bifactor)

pathLoadings = tmp$parameters[,c("row", "col", "Std.Estimate")]

pathLoadings$Std.Estimate = round(pathLoadings$Std.Estimate,2)

pathLoadings


                row              col Std.Estimate
1          Courtesy   general_factor         0.78
2      Friendliness   general_factor         0.76
3       Helpfulness   general_factor         0.84
4           Service   general_factor         0.81
5      Seat_Comfort   general_factor         0.79
6    Seat_Roominess   general_factor         0.72
7  Overhead_Storage   general_factor         0.70
8    Clean_Aircraft   general_factor         0.74
9  Easy_Reservation   general_factor         0.75
10  Preferred_Seats   general_factor         0.70
11   Flight_Options   general_factor         0.62
12    Ticket_Prices   general_factor         0.64
13 Easy_Reservation        ticketing         0.31
14  Preferred_Seats        ticketing         0.40
15   Flight_Options        ticketing         0.39
16    Ticket_Prices        ticketing         0.44
17     Seat_Comfort         aircraft         0.30
18   Seat_Roominess         aircraft         0.37
19 Overhead_Storage         aircraft         0.44
20   Clean_Aircraft         aircraft         0.34
21         Courtesy          service         0.27
22     Friendliness          service         0.44
23      Helpfulness          service         0.27
24          Service          service         0.34
25         Courtesy         Courtesy         0.33
26     Friendliness     Friendliness         0.23
27      Helpfulness      Helpfulness         0.23
28          Service          Service         0.23
29     Seat_Comfort     Seat_Comfort         0.28
30   Seat_Roominess   Seat_Roominess         0.34
31 Overhead_Storage Overhead_Storage         0.32
32   Clean_Aircraft   Clean_Aircraft         0.34
33 Easy_Reservation Easy_Reservation         0.34
34  Preferred_Seats  Preferred_Seats         0.36
35   Flight_Options   Flight_Options         0.47
36    Ticket_Prices    Ticket_Prices         0.39
> 

PS: If you want to learn more about fit indices, I recommend this Kenny link.

You probably want a figure as well, no? I describe doing that here.

R Code to Generate the Simulated Data and Run All Analyses

   
# First, we create a matrix of factor loadings.
# This pattern is called bifactor because it has a
# general factor with loadings from all the items
# and specific factors for separate components.
# The outcome variables are also formed as
# combinations of these general and specific factors.

loadings <- matrix(c(
 .33, .58, .00, .00,  # Ease of Making Reservation
 .35, .55, .00, .00,  # Availability of Preferred Seats
 .30, .52, .00, .00,  # Variety of Flight Options
 .40, .50, .00, .00,  # Ticket Prices
 .50, .00, .55, .00,  # Seat Comfort
 .41, .00, .51, .00,  # Roominess of Seat Area
 .45, .00, .57, .00,  # Availability of Overhead Storage
 .32, .00, .54, .00,  # Cleanliness of Aircraft
 .35, .00, .00, .50,  # Courtesy
 .38, .00, .00, .57,  # Friendliness
 .60, .00, .00, .50,  # Helpfulness
 .52, .00, .00, .58,  # Service
 .43, .10, .20, .30,  # Overall Satisfaction
 .35, .50, .40, .20,  # Purchase Intention
 .25, .50, .50, .00), # Willingness to Recommend
 nrow = 15, ncol = 4, byrow = TRUE
)

# Matrix multiplication produces the correlation matrix,
# except for the diagonal.
cor_matrix <- loadings %*% t(loadings)
# Diagonal set to ones.
diag(cor_matrix) <- 1
# cor_matrix = cov2cor(cor_matrix)
require(mvtnorm)
N = 1000
set.seed(7654321) #needed in order to reproduce the same data each time
std_ratings <- as.data.frame(rmvnorm(N, sigma = cor_matrix))

# Creates a mixture of two data sets:
# first 50 observations assigned uniformly lower scores.
ratings           <- data.frame(matrix(rep(0,15000),nrow=1000))
ratings[1:50,]    <- std_ratings[1:50,] * 2
ratings[51:1000,] <- std_ratings[51:1000,] * 2+7.0

# Ratings given different means
ratings[1] <- ratings[1]+2.2
ratings[2] <- ratings[2]+0.6
ratings[3] <- ratings[3]+0.3
ratings[4] <- ratings[4]+0.0
ratings[5] <- ratings[5]+1.5
ratings[6] <- ratings[6]+1.0
ratings[7] <- ratings[7]+0.5
ratings[8] <- ratings[8]+1.5
ratings[9] <- ratings[9]+2.4
ratings[10]<- ratings[10]+2.2
ratings[11]<- ratings[11]+2.1
ratings[12]<- ratings[12]+2.0
ratings[13]<- ratings[13]+1.5
ratings[14]<- ratings[14]+1.0
ratings[15]<- ratings[15]+0.5

# Truncates Scale to be between 1 and 9
ratings[ratings>9]<-9
ratings[ratings<1]<-1
# Rounds to single digit.
ratings<-round(ratings,0)

# Assigns names to the variables in the data frame called ratings
names(ratings) = c("Easy_Reservation","Preferred_Seats","Flight_Options","Ticket_Prices","Seat_Comfort","Seat_Roominess","Overhead_Storage","Clean_Aircraft","Courtesy","Friendliness","Helpfulness","Service","Satisfaction","Fly_Again","Recommend")