The most dependable way to read your SPSS file is to use Stat/Transfer to transform it to a Stata (.dta) file. In R, I use the readstata13 package to read Stata v13 files. I also loaded my statistics package (zStat). For convenience, I selected only the variables we need for these analyses. I added some variable labels for reference. Descriptive stats should be compared to SPSS to insure we are all on the same page.
library(readstata13)
library(zStat)
Brain1 = read.dta13('/prj/hnd/git/dataMgr/Projects/Beatty/2017-01-25-Brains/1/2017-04-12/Brain1-4.11.2017.dta')
Brain1 = Brain1[,zQ(LnLesion, Disc02b, TotDisc4, Age_3, SESComp, Sex)]
Brain1 = within(Brain1, {
zVarLbl(LnLesion) = 'ln(Lesion volume)'
zVarLbl(Disc02b) = 'Life harder'
zVarLbl(TotDisc4) = 'Racial discrimination'
})
zDesc(Brain1)
Variable N Mean SD Min Max Label
LnLesion LnLesion 71 6.04 1.63 0.0 9.68 ln(Lesion volume)
Disc02b Disc02b 71 2.00 1.01 1.0 4.00 Life harder
TotDisc4 TotDisc4 71 7.31 1.59 6.0 12.00 Racial discrimination
Age_3 Age_3 71 50.58 9.92 32.9 69.40
SESComp SESComp 71 0.55 0.50 0.0 1.00
Sex Sex 71 0.39 0.49 0.0 1.00
For complete sanity, re-run regression analyses to make sure all the coefficients agree. The regression results are stored in reg1 and reg2, which are R objects containing all the information we need to do the plots.
reg1 = lm(LnLesion ~ Disc02b*Age_3 + I(Disc02b^2)*Age_3 + SESComp + Sex, data=Brain1)
reg2 = lm(LnLesion ~ TotDisc4*Age_3 + I(TotDisc4^2)*Age_3 + SESComp + Sex, data=Brain1)
summary(reg1)
Call:
lm(formula = LnLesion ~ Disc02b * Age_3 + I(Disc02b^2) * Age_3 +
SESComp + Sex, data = Brain1)
Residuals:
Min 1Q Median 3Q Max
-5.5148969 -0.9449312 -0.0226640 0.9585741 3.1710691
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.40859070 4.84156405 -1.53021 0.1309718
Disc02b 11.64280855 5.10834238 2.27918 0.0260549
Age_3 0.27529426 0.09555082 2.88113 0.0054107
I(Disc02b^2) -2.54029234 1.09308191 -2.32397 0.0233658
SESComp 0.58582447 0.40928352 1.43134 0.1572746
Sex -0.19250979 0.39709106 -0.48480 0.6294996
Disc02b:Age_3 -0.24412679 0.09837437 -2.48161 0.0157588
Age_3:I(Disc02b^2) 0.05320830 0.02097933 2.53622 0.0136969
Residual standard error: 1.530283 on 63 degrees of freedom
Multiple R-squared: 0.2089847, Adjusted R-squared: 0.1210941
F-statistic: 2.377783 on 7 and 63 DF, p-value: 0.03190808
summary(reg2)
Call:
lm(formula = LnLesion ~ TotDisc4 * Age_3 + I(TotDisc4^2) * Age_3 +
SESComp + Sex, data = Brain1)
Residuals:
Min 1Q Median 3Q Max
-5.6727700 -0.6744797 -0.0485964 0.9986060 3.0518884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -96.705688577 29.806309661 -3.24447 0.00188495
TotDisc4 26.379859499 7.877274214 3.34886 0.00137369
Age_3 1.934106518 0.556942948 3.47272 0.00093679
I(TotDisc4^2) -1.669247822 0.501881927 -3.32598 0.00147306
SESComp 0.413772070 0.407984254 1.01419 0.31437309
Sex -0.362045021 0.383900932 -0.94307 0.34924826
TotDisc4:Age_3 -0.496357909 0.145991998 -3.39990 0.00117435
Age_3:I(TotDisc4^2) 0.031398524 0.009244291 3.39653 0.00118660
Residual standard error: 1.482345 on 63 degrees of freedom
Multiple R-squared: 0.2577672, Adjusted R-squared: 0.1752969
F-statistic: 3.125576 on 7 and 63 DF, p-value: 0.006762382
Find cut points for discrimination measures.
(mean1 = mean(Brain1$Disc02b))
[1] 2
(mean2 = mean(Brain1$TotDisc4))
[1] 7.309859155
(sd1 = sd(Brain1$Disc02b))
[1] 1.014185106
(sd2 = sd(Brain1$TotDisc4))
[1] 1.590970497
(Cuts1 = c(mean1-sd1, mean1, mean1+sd1))
[1] 0.9858148943 2.0000000000 3.0141851057
(Cuts2 = c(mean2-sd2, mean2, mean2+sd2))
[1] 5.718888658 7.309859155 8.900829652
Calculate predicted scores for plotting.
plot1 = zHat(Brain1, reg1, xAxis='Age_3', cutVars=list(Disc02b=Cuts1))
plot2 = zHat(Brain1, reg2, xAxis='Age_3', cutVars=list(TotDisc4=Cuts2))
zQuick(plot1)
Dimensions: 117 3
Age_3 Disc02b Hat
Min. :32 Min. :0.9858149 Min. :4.609077
1st Qu.:41 1st Qu.:0.9858149 1st Qu.:5.903299
Median :51 Median :2.0000000 Median :5.956053
Mean :51 Mean :2.0000000 Mean :6.074152
3rd Qu.:61 3rd Qu.:3.0141851 3rd Qu.:6.246272
Max. :70 Max. :3.0141851 Max. :7.889994
zQuick(plot2)
Dimensions: 117 3
Age_3 TotDisc4 Hat
Min. :32 Min. :5.718889 Min. :3.565237
1st Qu.:41 1st Qu.:5.718889 1st Qu.:6.013261
Median :51 Median :7.309859 Median :6.128923
Mean :51 Mean :7.309859 Mean :6.063655
3rd Qu.:61 3rd Qu.:8.900830 3rd Qu.:6.276972
Max. :70 Max. :8.900830 Max. :8.216482
# plot1
# plot2
Plot results for life harder
par(cex.axis=1.1, cex.lab=1.1, cex.main=1.25, las=1, lwd=2)
with(subset(plot1, Disc02b==Cuts1[1]), plot (Age_3, Hat, lty=2, type='l', xlim=c(30,70), xlab='Age', ylim=c(4,8), ylab='ln(Lesion volume)'))
with(subset(plot1, Disc02b==Cuts1[2]), lines(Age_3, Hat, lty=1))
with(subset(plot1, Disc02b==Cuts1[3]), lines(Age_3, Hat, lty=3))
legend(30, 8, zQ(High, Medium, Low), lty=c(3,1,2), bty='n')
Plot results for racial discrimination
par(cex.axis=1.1, cex.lab=1.1, cex.main=1.25, las=1, lwd=2)
with(subset(plot2, TotDisc4==Cuts2[1]), plot (Age_3, Hat, lty=2, type='l', xlim=c(30,70), xlab='Age', ylim=c(4,8), ylab='ln(Lesion volume)'))
with(subset(plot2, TotDisc4==Cuts2[2]), lines(Age_3, Hat, lty=1))
with(subset(plot2, TotDisc4==Cuts2[3]), lines(Age_3, Hat, lty=3))
legend(30, 8, zQ(High, Medium, Low), lty=c(3,1,2), bty='n')
These plots might be more interesting if discrimination were on the x-axis with age cuts.
Cuts = c(40,50,60)
plot1 = zHat(Brain1, reg1, xAxis='Disc02b', cutVars=list(Age_3=Cuts), xIncr=.1)
plot2 = zHat(Brain1, reg2, xAxis='TotDisc4', cutVars=list(Age_3=Cuts), xIncr=.1)
zQuick(plot1)
Dimensions: 93 3
Disc02b Age_3 Hat
Min. :1.0 Min. :40 Min. :4.768632
1st Qu.:1.7 1st Qu.:40 1st Qu.:5.900123
Median :2.5 Median :50 Median :5.977985
Mean :2.5 Mean :50 Mean :6.040031
3rd Qu.:3.3 3rd Qu.:60 3rd Qu.:6.131653
Max. :4.0 Max. :60 Max. :7.771029
zQuick(plot2)
Dimensions: 183 3
TotDisc4 Age_3 Hat
Min. : 6.0 Min. :40 Min. :-0.4665957
1st Qu.: 7.5 1st Qu.:40 1st Qu.: 5.5489703
Median : 9.0 Median :50 Median : 6.0731899
Mean : 9.0 Mean :50 Mean : 5.7888686
3rd Qu.:10.5 3rd Qu.:60 3rd Qu.: 6.3515951
Max. :12.0 Max. :60 Max. : 9.5173844
# plot1
# plot2
Plot results for life harder
par(cex.axis=1.1, cex.lab=1.1, cex.main=1.25, las=1, lwd=2)
with(subset(plot1, Age_3==Cuts[1]), plot (Disc02b, Hat, lty=2, type='l', xlim=c(1,4), ylim=c(4,10), xlab='Life harder', ylab='ln(Lesion volume)'))
with(subset(plot1, Age_3==Cuts[2]), lines(Disc02b, Hat, lty=1))
with(subset(plot1, Age_3==Cuts[3]), lines(Disc02b, Hat, lty=3))
legend(1, 10, zQ('40yo', '50yo', '60yo'), lty=c(2,1,3), bty='n')
Plot results for racial discrimination
par(cex.axis=1.1, cex.lab=1.1, cex.main=1.25, las=1, lwd=2)
with(subset(plot2, Age_3==Cuts[1]), plot (TotDisc4, Hat, lty=2, type='l', xlim=c(6,12), ylim=c(0,10), xlab='Racial discrimination', ylab='ln(Lesion volume)'))
with(subset(plot2, Age_3==Cuts[2]), lines(TotDisc4, Hat, lty=1))
with(subset(plot2, Age_3==Cuts[3]), lines(TotDisc4, Hat, lty=3))
legend(6, 10, zQ('40yo', '50yo', '60yo'), lty=c(2,1,3), bty='n')
/prj/hnd/git/dataMgr/Projects/Beatty/2017-01-25-Brains/1/2017-04-12/2017-04-12.Rmd