The naive approach would assume all the responses from all subjects are independent and estimate the parameters and standard errors with the usual logistic model, which can be done with R's glm function.

library(geepack)
gss = read.csv(file.choose(), header=T) # "gss.confidence.csv"
# all observations independent
fit.glm = glm(greatly~question, data=gss, family=binomial(link='logit'))

By default, R will treat the question about education (first alphabetically among the question types) as the baseline. Compare the output below to the estimates and standard errors obtained from the calculations above.

> summary(fit.glm)

Call:
glm(formula = greatly ~ question, family = binomial(link = "logit"), 
    data = gss)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.099  -0.941  -0.762   1.258   1.660  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -1.0868     0.0601  -18.08   <2e-16 ***
questionconmedic   0.5022     0.0811    6.19    6e-10 ***
questionconsci     0.8995     0.0798   11.27   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5722.8  on 4400  degrees of freedom
Residual deviance: 5591.4  on 4398  degrees of freedom
AIC: 5597

If instead, we estimate the parameters and standard errors according to the GEE approach, we use the geeglm function and have several choices for the working correlation structure (within-subjects). Note that even with the working correlation structure specified as independent, as it is below with the corstr option, the sample correlations among the three responses within subjects are still being incorporated into the standard error estimates. 

# gee with independent corr structure 
fit.ind = geeglm(greatly~question, data=gss, id=id, family=binomial, corstr="independence", scale.fix=T)

The id= option is where we specify the clusters within which we have repeated observations that may be correlated, and tehe scale.fix=T option is to avoid R's default approach of introducing an overdispersion scale parameter. Compared with the glm function output above, the estimates are identical, but the standard errors are slightly smaller. The ANOVA function can also be used to output an ANOVA table for an overall effect due to question type. Versus the reference chi-square distribution with two degrees of freedom, the test statistic of 149 is strong evidence that subjects' confidence differs among the three institutions, with the highest proportion of "A great deal" of confidence applying to the scientific community.

> summary(fit.ind)

Call:
geeglm(formula = greatly ~ question, family = binomial, data = gss, 
    id = id, corstr = "independence", scale.fix = TRUE)

 Coefficients:
                 Estimate Std.err Wald Pr(>|W|)    
(Intercept)       -1.0868  0.0601  327  < 2e-16 ***
questionconmedic   0.5022  0.0697   52  5.6e-13 ***
questionconsci     0.8995  0.0737  149  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = independence 
Scale is fixed.

Number of clusters:   1467  Maximum cluster size: 3 

> anova(fit.ind)
Analysis of 'Wald statistic' Table
Model: binomial, link: logit
Response: greatly
Terms added sequentially (first to last)

         Df  X2 P(>|Chi|)    
question  2 149    <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Other options for the working correlation structure can produce slightly different results, but in general, the results are not that sensitive to the specified working correlation structure because empirical correlations among the responses from the data are largely dominant in the GEE calculations. With this in mind, a good compromise between model fit and parameter sparsity is the exchangeable structure, which allows a single correlation parameter for all pairwise responses within a subject. The results below are for the exchangeable structure; these are particularly similar to those assuming the independent correlation structure in this example because the estimated correlation in the data is so small.

# gee with exchangeable corr structure
fit.exch = geeglm(greatly~question, data=gss, id=id, family=binomial, corstr="exchangeable", scale.fix=T)
# gee with unstructured corr structure
fit.un = geeglm(greatly~question, data=gss, id=id, family=binomial, corstr="unstructured", scale.fix=T)

> summary(fit.exch)

Call:
geeglm(formula = greatly ~ question, family = binomial, data = gss, 
    id = id, corstr = "exchangeable", scale.fix = TRUE)

 Coefficients:
                 Estimate Std.err Wald Pr(>|W|)    
(Intercept)       -1.0868  0.0601  327  < 2e-16 ***
questionconmedic   0.5022  0.0697   52  5.6e-13 ***
questionconsci     0.8995  0.0737  149  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = exchangeable 
Scale is fixed.

  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha    0.235  0.0173
Number of clusters:   1467  Maximum cluster size: 3  

# gee with age covariate
gss$age = gss$age - mean(gss$age)
fit.age = geeglm(greatly~question+age, data=gss, id=id, family=binomial, corstr="exchangeable", scale.fix=T)

We see that, in addition to the institution in question, a person's confidence is also related with age. For each additional year in age, the estimated odds of having "A great deal" of confidence in the education system (or either of the others) is multiplied by

\( e^{-0.00504}=0.995\)

So, the proportions of people having "A great deal" of confidence in these institutions tend to be smaller for older ages, but this effect is very small in practical terms. The statistical significance is likely due to the large sample size.

> summary(fit.age)

Call:
geeglm(formula = greatly ~ question + age, family = binomial, 
    data = gss, id = id, corstr = "exchangeable", scale.fix = T)

 Coefficients:
                 Estimate  Std.err   Wald Pr(>|W|)    
(Intercept)      -1.08903  0.06026 326.65  < 2e-16 ***
questionconmedic  0.50306  0.06979  51.95  5.7e-13 ***
questionconsci    0.90109  0.07381 149.06  < 2e-16 ***
age              -0.00504  0.00228   4.91    0.027 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = exchangeable 
Scale is fixed.

  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha    0.233  0.0174
Number of clusters:   1467  Maximum cluster size: 3