Next: About this document
Up: No Title
Previous: Data program example:
R : Copyright 2001, The R Development Core Team
Version 1.3.1 (2001-08-31)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.
R is a collaborative project with many contributors.
Type `contributors()' for more information.
Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.
> invisible(options(echo = TRUE))
> # R
> # using data from year 1996
>
> # CF0301 is 7-PT SCALE PARTY IDENTIFICATION
> # CF0704 is PARTY OF PRES VOTE- ALL MAJOR CANDIDATES
> # CF0704A is PARTY OF PRES VOTE- 2 MAJOR PARTIES
> # CF0724 is DID R WATCH TV PROGRAMS ABOUT CAMPAIGN
> # VCF0803 is LIBERAL-CONSERVATIVE 7PT SCALE
> # CF9088 is DEM PRES CAND- 7PT LIBERAL/CONSERV SCALE
> # CF9096 is REP PRES CAND- 7PT LIBERAL/CONSERV SCALE
> DF <- read.table(file="tmp.sas2awk.dat",header=T);
> attach(DF);
>
> CF0301 <- ifelse ( CF0301 == 0 , NA , CF0301 );
> CF0704 <- ifelse ( CF0704 == 0 , NA , CF0704 );
> CF0704A <- ifelse ( CF0704A == 0 , NA , CF0704A );
> CF0704A <- ifelse ( CF0704A == 0 , NA , CF0704A );
> CF0724 <- ifelse ( CF0724 == 0 , NA , CF0724 );
> VCF0803 <- ifelse ( VCF0803 == 0 , NA , VCF0803 );
> CF9088 <- ifelse ( CF9088 == 0 , NA , CF9088 );
> CF9088 <- ifelse ( CF9088 >= 8 , NA , CF9088 );
> CF9096 <- ifelse ( CF9096 == 0 , NA , CF9096 );
> CF9096 <- ifelse ( CF9096 >= 8 , NA , CF9096 );
>
> # get rid of unwanted value in party ID
> pid <- CF0301;
> pid <- ifelse(pid == 9, NA, pid);
>
> # get rid of unwanted value in libcon self-placement
> libcon <- VCF0803;
> libcon <- ifelse(libcon == 9, NA, libcon);
>
> # presvote3 is three major party candidates only
> presvote3 <- CF0704;
>
> # presvote is Dem and Rep candidates only
> presvote <- CF0704A - 1;
>
> # watchtv is how much watched campaign on TV
> watchtv <- CF0724;
>
> # compute spatial model measure of difference between Dem and Rep candidates
> demlibcon <- CF9088;
> replibcon <- CF9096;
> demdist <- abs(libcon-demlibcon);
> repdist <- abs(libcon-replibcon);
> distdiff <- demdist-repdist;
>
> # crosstab of libcon given pid
> XtableXXX <- xtabs(wgtXXX ~ libcon + pid , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
pid
libcon 1 2 3 4 5 6 7
1 11.0338 6.0570 4.8070 1.0021 0.0000 0.0000 0.0000
2 65.2037 29.2879 19.0274 2.9781 3.7373 3.1289 1.9285
3 46.4409 61.2765 42.2474 8.4403 5.1316 14.0725 1.9081
4 68.2672 88.9166 76.8295 50.0696 43.5462 55.0045 21.9070
5 22.4832 38.1181 27.1080 19.2086 49.1673 71.2493 30.4313
6 19.3470 18.6181 7.5078 6.1409 36.5094 61.7322 127.2785
7 2.9321 2.1205 2.0426 1.2132 7.0613 7.3538 21.0486
> # independence test using data from year 1996
> summary(XtableXXX);
Call: xtabs(formula = wgtXXX ~ libcon + pid, drop.unused.levels = T)
Number of cases in table: 1310.921
Number of factors: 2
Test for independence of all factors:
Chisq = 631.4, df = 36, p-value = 6.976e-110
Chi-squared approximation may be incorrect
> # weighted column proportions using data from year 1996
> round(t(t(XtableXXX) / apply(XtableXXX,2,sum)), digits=3);
pid
libcon 1 2 3 4 5 6 7
1 0.047 0.025 0.027 0.011 0.000 0.000 0.000
2 0.277 0.120 0.106 0.033 0.026 0.015 0.009
3 0.197 0.251 0.235 0.095 0.035 0.066 0.009
4 0.290 0.364 0.428 0.562 0.300 0.259 0.107
5 0.095 0.156 0.151 0.216 0.339 0.335 0.149
6 0.082 0.076 0.042 0.069 0.252 0.290 0.622
7 0.012 0.009 0.011 0.014 0.049 0.035 0.103
>
> # crosstab of libcon given pid and watchtv
> XtableXXX <- xtabs(wgtXXX ~ libcon + pid + watchtv , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
, , watchtv = 1
pid
libcon 1 2 3 4 5 6 7
1 2.9225 2.1635 0.0000 0.0000 0.0000 0.0000 0.0000
2 12.5496 10.8378 7.9023 1.9781 0.8794 0.9577 0.9081
3 3.9002 10.9099 9.5101 2.9735 0.9899 4.1070 0.0000
4 8.7117 18.4375 17.2989 15.7230 13.6720 11.1897 4.1729
5 1.0000 11.8085 4.1371 6.5478 8.1683 13.2768 4.8868
6 2.2334 3.0630 2.7518 0.9081 4.2074 11.1635 13.6575
7 0.0000 1.1205 1.0000 1.2132 1.2902 0.0000 0.9081
, , watchtv = 2
pid
libcon 1 2 3 4 5 6 7
1 6.1529 2.8935 3.8070 1.0021 0.0000 0.0000 0.0000
2 50.5943 16.4998 8.3944 1.0000 2.8579 0.0000 1.0204
3 36.9960 42.6713 24.5909 5.4668 2.9285 7.7966 1.9081
4 48.6916 62.8219 53.5414 30.1307 23.8291 39.6604 15.8027
5 20.4832 22.2418 18.0512 11.4915 37.0967 51.4836 24.6372
6 16.2241 14.3146 3.7560 5.2328 28.6872 45.6148 107.7144
7 2.9321 1.0000 0.0000 0.0000 3.6167 6.1826 16.6216
> # weighted column proportions using data from year 1996
> XnXXX <- length(dim(XtableXXX));
> XtsumXXX <- apply(XtableXXX,2:XnXXX,sum);
> XaXXX <- array(XtsumXXX,dim=dim(XtableXXX)[c(2:XnXXX,1)]);
> round(XtableXXX/aperm(XaXXX,c(XnXXX,1:(XnXXX-1))), digits=3);
, , watchtv = 1
pid
libcon 1 2 3 4 5 6 7
1 0.093 0.037 0.000 0.000 0.000 0.000 0.000
2 0.401 0.186 0.185 0.067 0.030 0.024 0.037
3 0.125 0.187 0.223 0.101 0.034 0.101 0.000
4 0.278 0.316 0.406 0.536 0.468 0.275 0.170
5 0.032 0.202 0.097 0.223 0.280 0.326 0.199
6 0.071 0.053 0.065 0.031 0.144 0.274 0.557
7 0.000 0.019 0.023 0.041 0.044 0.000 0.037
, , watchtv = 2
pid
libcon 1 2 3 4 5 6 7
1 0.034 0.018 0.034 0.018 0.000 0.000 0.000
2 0.278 0.102 0.075 0.018 0.029 0.000 0.006
3 0.203 0.263 0.219 0.101 0.030 0.052 0.011
4 0.267 0.387 0.477 0.555 0.241 0.263 0.094
5 0.112 0.137 0.161 0.212 0.375 0.342 0.147
6 0.089 0.088 0.033 0.096 0.290 0.303 0.642
7 0.016 0.006 0.000 0.000 0.037 0.041 0.099
>
> # crosstab of three-party vote given pid and watchtv
> XtableXXX <- xtabs(wgtXXX ~ presvote3 + pid , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
pid
presvote3 1 2 3 4 5 6 7
1 237.8780 166.7259 97.9056 21.8474 23.0423 34.8925 9.3169
2 4.9358 17.2297 7.5622 22.3237 77.5199 119.9529 176.7608
3 6.1702 13.4894 20.4324 10.3412 12.7375 17.2742 1.9285
> # independence test using data from year 1996
> summary(XtableXXX);
Call: xtabs(formula = wgtXXX ~ presvote3 + pid, drop.unused.levels = T)
Number of cases in table: 1100.267
Number of factors: 2
Test for independence of all factors:
Chisq = 702.7, df = 12, p-value = 1.137e-142
Chi-squared approximation may be incorrect
> # weighted column proportions using data from year 1996
> round(t(t(XtableXXX) / apply(XtableXXX,2,sum)), digits=3);
pid
presvote3 1 2 3 4 5 6 7
1 0.955 0.844 0.778 0.401 0.203 0.203 0.05
2 0.020 0.087 0.060 0.410 0.684 0.697 0.94
3 0.025 0.068 0.162 0.190 0.112 0.100 0.01
>
> # crosstab of two-party vote given pid and watchtv
> XtableXXX <- xtabs(wgtXXX ~ presvote + pid , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
pid
presvote 1 2 3 4 5 6 7
0 237.8780 166.7259 97.9056 21.8474 23.0423 34.8925 9.3169
1 4.9358 17.2297 7.5622 22.3237 77.5199 119.9529 176.7608
> # independence test using data from year 1996
> summary(XtableXXX);
Call: xtabs(formula = wgtXXX ~ presvote + pid, drop.unused.levels = T)
Number of cases in table: 1017.894
Number of factors: 2
Test for independence of all factors:
Chisq = 639.3, df = 6, p-value = 7.805e-135
> # weighted column proportions using data from year 1996
> round(t(t(XtableXXX) / apply(XtableXXX,2,sum)), digits=3);
pid
presvote 1 2 3 4 5 6 7
0 0.98 0.906 0.928 0.495 0.229 0.225 0.05
1 0.02 0.094 0.072 0.505 0.771 0.775 0.95
>
> # one-way frequency table of libcon placements
> XtableXXX <- xtabs(wgtXXX ~ libcon , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
libcon
1 2 3 4 5 6 7
23.8465 126.3172 180.5173 406.6769 257.7658 277.1339 45.7721
> # weighted proportions using data from year 1996
> round(XtableXXX / sum(XtableXXX), digits=3);
libcon
1 2 3 4 5 6 7
0.018 0.096 0.137 0.309 0.196 0.210 0.035
> XtableXXX <- xtabs(wgtXXX ~ demlibcon , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
demlibcon
1 2 3 4 5 6 7
170.3448 468.4437 353.7143 340.6165 124.0586 96.7575 48.6063
> # weighted proportions using data from year 1996
> round(XtableXXX / sum(XtableXXX), digits=3);
demlibcon
1 2 3 4 5 6 7
0.106 0.292 0.221 0.213 0.077 0.060 0.030
> XtableXXX <- xtabs(wgtXXX ~ replibcon , drop.unused.levels = T);
> # weighted frequencies using data from year 1996
> XtableXXX;
replibcon
1 2 3 4 5 6 7
38.3328 61.5420 107.4572 222.2778 318.7721 634.9199 171.0126
> # weighted proportions using data from year 1996
> round(XtableXXX / sum(XtableXXX), digits=3);
replibcon
1 2 3 4 5 6 7
0.025 0.040 0.069 0.143 0.205 0.408 0.110
>
> # simple ordinary least squares regression model of two-party vote
> # explanatory variables are PID dummy variables and spatial model difference
> summary(lm(presvote ~ factor(pid) + distdiff, weights=wgtXXX));
Call:
lm(formula = presvote ~ factor(pid) + distdiff, weights = wgtXXX)
Residuals:
Min 1Q Median 3Q Max
-0.89447 -0.08392 0.01389 0.12353 1.00832
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.162203 0.023292 6.964 6.78e-12 ***
factor(pid)2 0.037511 0.031154 1.204 0.229
factor(pid)3 -0.010208 0.036115 -0.283 0.778
factor(pid)4 0.392531 0.056032 7.006 5.12e-12 ***
factor(pid)5 0.490806 0.040472 12.127 < 2e-16 ***
factor(pid)6 0.534850 0.037271 14.350 < 2e-16 ***
factor(pid)7 0.562966 0.039946 14.093 < 2e-16 ***
distdiff 0.071434 0.005446 13.116 < 2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.269 on 823 degrees of freedom
Multiple R-Squared: 0.7061, Adjusted R-squared: 0.7036
F-statistic: 282.5 on 7 and 823 DF, p-value: 0
>
> # simple probit regression model of two-party vote
> # explanatory variables are PID dummy variables and spatial model difference
> summary(glm(presvote ~ factor(pid) + distdiff,
+ family=binomial(link="probit"), weights=wgtXXX));
Call:
glm(formula = presvote ~ factor(pid) + distdiff, family = binomial(link = "probit"),
weights = wgtXXX)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.51471 -0.22324 -0.02203 0.25389 3.32229
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.60873 0.24270 -6.629 3.39e-11 ***
factor(pid)2 0.58885 0.29706 1.982 0.0474 *
factor(pid)3 0.23809 0.35097 0.678 0.4975
factor(pid)4 1.61944 0.35836 4.519 6.21e-06 ***
factor(pid)5 1.79791 0.29994 5.994 2.04e-09 ***
factor(pid)6 2.02980 0.28556 7.108 1.18e-12 ***
factor(pid)7 2.31723 0.31870 7.271 3.57e-13 ***
distdiff 0.46999 0.04687 10.027 < 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: 1124.56 on 830 degrees of freedom
Residual deviance: 360.83 on 823 degrees of freedom
AIC: 388.56
Number of Fisher Scoring iterations: 5
Warning message:
non-integer #successes in a binomial glm! in: eval(expr, envir, enclos)
> proc.time()
[1] 5.86 0.57 6.36 0.00 0.00
>
Walter Mebane
Tue Nov 13 18:27:23 EST 2001