Download Analysis of Learning Curves and Proportional Odds Model in Educational Psychology and more Assignments Biostatistics in PDF only on Docsity!
Biostat/Stat 571 Exercise
Answer Key
January, 2007
Question 2:
TenHave and Uttal (1994) report on the results of a psychology experiment in which children were shown a map giving the location of a hidden toy, and allowed up to three chances to successfully find the object. The goal of the experiment is to compare a group of children presented with a correctly oriented map (control) to a second group that was shown a rotated map. Are the children presented with the rotated map able to compensate and find the hidden object as well as the control children? Each child repeated the test 10 times. During each test, a child will only attempt the maze a second time if they fail the first time. Similarly, a child will only attempt the maze a third time if they fail the second time. If a child is unsuccessful the third time the outcome is an overall failure. Thus the outcome, number of attempts to successfully find the object takes (Y ), takes on one of four values with Y = 4 denoting all three attempts having failed.
(a) Figure (1) displays the distribution of the outcome, at each of the ten tests, by treatment group.
1 2 3 4 5 6 7 8 9 10
0
20
40
60
80
100
Test number
Percentage
Fail Third Second First
(a) Control Group (n = 46)
1 2 3 4 5 6 7 8 9 10
0
20
40
60
80
100
Test number
Percentage
Fail Third Second First
(b) Rotated Group (n = 43)
Figure 1: Percentage of the number of attempts per test
It is difficult to discern a “learning effect” over time for the group that were presented the correctly oriented map (control). It does not seem as though the percentage that fail decreases over time, nor that the percentage that pass 1st^ time increases with time. In the treatment group, that were presented with the rotated maps, it seems at though there is some evidence of a “learning effect”. In particular, the percentage that fail decreases across the 10 tests, while the percentage that pass 1st time increases.
(b) Let ∆i = Yi 8 − Yi 1 denote the difference between the outcomes of the 8th^ and 1st^ tests for the ith child. Consequently, a negative change indicates an improvement over time. For example, ∆ = − 3
corresponds to a score of 1 on the 8th^ test and a score of 4 on the 1st, indicating improvement. The table below summarises these differences by treatment group.
∆ -3 -2 -1 0 1 2 3 Total Rotated 10 14 5 8 3 2 1 43 Control 5 2 13 18 5 1 2 46
From this we find that 63% (29/43) of the rotated group improved their score from the 1st^ to the 8 th^ tests, with 10 children passing on the first attempt during the 8th^ test when they had failed all three attempts during the 1st^ test. In the control group, only 47% improved their score with a much larger group (42%) maintaining their score between the two tests. It is interesting to note that of the 8 children in the rotate group that maintained their score, 7 maintained a score of 4. However, in the control group of the 18 that maintained their score 16 maintained scores of 1 or 2. Neither treatment group exhibited much learning between the two tests. Due to the sparse data in the final two outcome categories, these have been collapsed into one category (although the results do not change considerably keeping the two seperate).
Score test for the proportional odds assumption (provided by SAS) yields a test statistic of 14.9 on 4 degrees of freedom (p-value = 0.0048), which suggests evidence against the null hypothesis that the assumption holds. Thus we should consider fitting seperate logisitc regressions for each dichotomisation of the outcome. These are presented (along with the fit of the proportional odds model) below.
I[∆i≤−3] I[∆i≤−2] I[∆i≤−1] I[∆i≤0] I[∆i ≤1] Ordinal(a)
Intercept(−3) -2.104 -2. Intercept(−2) -1.717 -1. Intercept(−1) -0.262 -0. Intercept(0) 1.588 1. Intercept(1) 2.662 2.
Rotate 0.910 1.951 0.990 0.261 -0.072 1.
(a) (^) POM: logit(γij ) = θj + Rotatei βR
In this model, for this formulation, we interpret the ‘Rotate’ coefficient in terms of the odds of having a “low” change (which would suggest learning) versus having a “high” change (which would suggest getting worse). For the logisitc model with outcome I[∆i≤−1], the odds of a learning effect versus no learning effect/worsening are estimated to be 2.69 times higher (exp^0.^99 ≈ 2.69; 95% CI (1.14, 6.39)) in the rotated group compared to the control group. This suggests that the rotated group is “more likely” to experience a learning effect (i.e. improvement in score).
(c) To examine the change between Test 1 and Test 8, we could also consider modelling the score at for Test 8 as a function of the score at time 1 as well as the treatment group. Initially, an interation between the Test 1 score and the treatment group was considered. The resulting likelihood ratio test for the interaction terms, in the proportional odds model, yields a test statistic of 7.23 on 3 degrees of freedom (p-value = 0.065). Further examination of the cofficients indicated that this was completely driven by one parameter and consequently the interactions were not incorporated. Note, this parameter was associated with the group that had a Test 1 score of 3, which had the lowest cell
Control Rotated
Test 8 score Test 8 Score 1 2 3 4 1 2 3 4 1 13 2 1 2 1 1 0 0 1 Test 1 2 9 3 2 0 2 2 0 2 2 Score 3 1 4 1 1 3 4 0 0 1 4 5 1 0 1 4 10 10 3 7
Consequently, we can see that trying to characterise/assess any learning effect is difficult (at best!) with this formulation.
S-Plus Code:
##########################
Exercise 1.
- Question 2.
##########################
Preliminary set-up; including loadin the MASS functions and Design library
remove( ls() ) options( contrasts = c("contr.treatment", "contr.poly") ) source( "polr.q" ) source( "vcov.q" ) source( "misc.q" ) library( Design )
###################################################################################
Load in the tenhave data, label variables and make all data manipulations
###################################################################################
tenhave <- matrix( scan( "tenhave.data", sep = " " ), ncol = 11, byrow = T ) dimnames( tenhave ) <- list( NULL, c( "rotate", paste( "test", 1:10, sep = "" ) ) ) tenhave <- as.data.frame( tenhave )
Change the coding of the ’rotate’ variable so that rotate = 1 is the rotated (intervention) group
tenhave$rotate <- 1 - tenhave$rotate
Compute ‘change’ outcome; Test8 - Test
- positive change implies improvement over time
tenhave$change <- tenhave$test8 - tenhave$test
##################
Part (a)
##################
Look at some individual ’learning curves’
sample.rotate <- sort( sample( c(1:43), 8 ) ) sample.not <- sort( sample( c(44:89), 8 ) ) par( mfrow = c(4,4) ) for( i in 1:8 ) { plot( 1:10, tenhave[sample.rotate[i], 2:11], type = "l", ylim = c(0, 4), xlab = "Test number", ylab = "Outcome", main = paste("ID:", sample.rotate[i], "- rotated") ) } for( i in 1:8 ) {
plot( 1:10, tenhave[sample.not[i], 2:11], type = "l", ylim = c(0, 4), xlab = "Test number", ylab = "Outcome", main = paste("ID:", sample.not[i], "- not") ) }
Calculate percentages within each test; PBT = percent by test
pbt.rotated <- matrix( 0, ncol = 10, nrow = 4 ) # n = 43 pbt.not <- matrix( 0, ncol = 10, nrow = 4 ) # n = 46 for( i in 1:10 ) { pbt.rotated[,i] <- (table( tenhave[1:43,(i+1)] ) / 43) * 100 pbt.not[,i] <- (table( tenhave[44:89,(i+1)] ) / 46) * 100 }
Produce barplots
barplot( pbt.not, names = as.character(c(1:10)), xlab = "Test number", ylab = "Percentage", col = c(1,5,7,6) ) legend( 8.7, 24, c("Fail", "Third", "Second", "First" ), fill = c(6,7,5,1) ) barplot( pbt.rotated, names = as.character(c(1:10)), xlab = "Test number", ylab = "Percentage", col = c(1,5,7,6) ) legend( 0.5, 97, c("Fail", "Third", "Second", "First" ), fill = c(6,7,5,1) )
##################
Part (b)
##################
table( tenhave$rotate, tenhave$change ) table( tenhave$test1[tenhave$rotate == 0 & tenhave$change == 0], tenhave$test8[tenhave$rotate == 0 & tenhave$change == 0] ) table( tenhave$test1[tenhave$rotate == 1 & tenhave$change == 0], tenhave$test8[tenhave$rotate == 1 & tenhave$change == 0] )
Sparse outcome for last two categories; collapse
- should be ok since under the PO assumption regression parameter is the same
- only the intercepts will change
tenhave$change[tenhave$change == 3] <- 2
Perform a series of logistic regressions by sequentially changing the dichotomisation
- out1 : I_{change \le -3}
- out2 : I_{change \le -2}
.....
- out5 : I_{change \le 1}
- out6 : I_{change \le 2} ## not used
tenhave$out1 <- ifelse( (tenhave$change < -2), 1, 0 ) tenhave$out2 <- ifelse( (tenhave$change < -1), 1, 0 ) tenhave$out3 <- ifelse( (tenhave$change < 0), 1, 0 ) tenhave$out4 <- ifelse( (tenhave$change < 1), 1, 0 ) tenhave$out5 <- ifelse( (tenhave$change < 2), 1, 0 )
tenhave$out6 <- ifelse( (tenhave$change < 3), 1, 0 )
model.out1 <- glm( out1 ~ rotate, data = tenhave, family = binomial ) model.out2 <- glm( out2 ~ rotate, data = tenhave, family = binomial ) model.out3 <- glm( out3 ~ rotate, data = tenhave, family = binomial ) model.out4 <- glm( out4 ~ rotate, data = tenhave, family = binomial ) model.out5 <- glm( out5 ~ rotate, data = tenhave, family = binomial )
model.out6 <- glm( out6 ~ rotate, data = tenhave, family = binomial )
coef.output.b <- c( summary(model.out1)$coef[1,1:2], summary(model.out1)$coef[2,1:2] ) coef.output.b <- rbind( coef.output.b, c( summary(model.out2)$coef[1,1:2], summary(model.out2)$coef[2,1:2] ) ) coef.output.b <- rbind( coef.output.b, c( summary(model.out3)$coef[1,1:2], summary(model.out3)$coef[2,1:2] ) ) coef.output.b <- rbind( coef.output.b, c( summary(model.out4)$coef[1,1:2], summary(model.out4)$coef[2,1:2] ) ) coef.output.b <- rbind( coef.output.b, c( summary(model.out5)$coef[1,1:2], summary(model.out5)$coef[2,1:2] ) )
coef.output.b <- rbind( coef.output.b, c( summary(model.out6)$coef[1,1:2], summary(model.out6)$coef[2,1:2] ) )
coef.output.b <- as.data.frame( round( coef.output.b, 3 ) ) names( coef.output.b ) <- c( "Intercept", "se", "rotate", "se" ) coef.output.b
Proportional Odds Model for change
model.POMb <- polr( factor(change) ~ rotate, data = tenhave )
elements of Σ. This is because Cov(B(Y − π)) = BΣBT^ : the negative off-diagonals of Σ (a multinomial covariance) are multiplied by the positive off-diagonals of B, and subtracted from Var(Y), since the diagonal elements of B are 1. 2c) Using code adapted from the WESDR example (page 67 of the class notes):
> tenhave[1:3,] rotate t1 t2 t3 t4 t5 t6 t7 t8 t9 t 1 1 4 4 4 4 4 4 4 2 4 4 2 1 2 4 4 4 2 1 2 4 3 4 3 1 4 4 4 2 4 4 3 2 4 4
construct continuation model data for trials 1 and 8
(1) construct pairs (Y,H)
nsubj_ nrow(tenhave) ncuts_ 3 y1.crm _NULL y8.crm _NULL h1.crm _NULL h8.crm _NULL id _NULL for (j in (1:nsubj)) { yj1 _rep(0,ncuts) yj8 _rep(0,ncuts) if (tenhave$t1[j] <= ncuts ) yj1[tenhave$t1[j]] _ 1 if (tenhave$t8[j] <= ncuts ) yj8[tenhave$t8[j]] _ 1 hj1 _ 1-c(0,cumsum(yj1)[1:(ncuts-1)]) hj8 _ 1-c(0,cumsum(yj8)[1:(ncuts-1)]) y1.crm _ c(y1.crm, yj1) y8.crm _ c(y8.crm, yj8) h1.crm _ c(h1.crm, hj1) h8.crm _ c(h8.crm, hj8) id _ c(id, rep(j, ncuts)) }
(2) construct intercepts
level _ factor(rep(1:ncuts,nsubj))int.mat _ NULL for (j in 1:ncuts) { intj _ rep(0,ncuts) intj[j] _ 1 int.mat _ cbind(int.mat, rep(intj, nsubj)) }
(3) expand the covariate
rotate.crm _ rep(tenhave$rotate, rep(ncuts, nsubj))
(4) drop H=0 and build dataframe
keep1 _ h1.crm== keep8 _ h8.crm== t1.crm.data _ data.frame( id=id[keep1], y=y1.crm[keep1], level=level[keep1], rotate=rotate.crm[keep1]) t8.crm.data _ data.frame( id=id[keep8], y=y8.crm[keep8], level=level[keep8], rotate=rotate.crm[keep8])
model fitting for 2c
fit1 _ glm(y~level+rotate, family=binomial, data=t1.crm.data) summary(fit1) Coefficients: Value Std. Error t value (Intercept) -0.5442677 0.2831311 -1. level2 0.6443691 0.4058320 1. level3 0.6612261 0.4712200 1. rotate -2.0075675 0.3800151 -5.
fit2 _ glm(y~level*rotate, family=binomial, data=t1.crm.data) summary(fit2) Coefficients: Value Std. Error t value (Intercept) -0.4418328 0.3021090 -1. level2 0.4418328 0.4838667 0. level3 0.4418328 0.6139903 0. rotate -2.5763281 0.7634178 -3. level2rotate 0.8127395 0.9596033 0. level3rotate 0.7845687 1.0496489 0.
anova(fit1,fit2)Analysis of Deviance Table
Response: y Terms Resid. Df Resid. Dev Test Df Deviance 1 level + rotate 203 199. 2 level * rotate 201 198.8258 +level:rotate 2 0.
An analysis of deviance indicates that a common treatment odds ratio is appropriate (1-pchisq(.859,2)= 0.65).
model fitting for 2d
fit1d _ glm(y~level+rotate, family=binomial, data=t8.crm.data) summary(fit1d) Coefficients: Value Std. Error t value (Intercept) 0.4130039 0.2688397 1. level2 -0.1236940 0.3780036 -0. level3 -0.4010700 0.4817271 -0. rotate -0.8071621 0.3301882 -2.
fit2d _ glm(y~level*rotate, family=binomial, data=t8.crm.data) summary(fit2d) Coefficients: Value Std. Error t value (Intercept) 0.44182741 0.3018794 1. level2 -0.21868390 0.5622350 -0. level3 -0.44182741 0.7688505 -0. rotate -0.86670652 0.4339232 -1. level2rotate 0.17356755 0.7585855 0. level3rotate 0.07851702 0.9874841 0.
anova(fit1d, fit2d) Analysis of Deviance Table Response: y
1 level + rotateTerms Resid. Df Resid. Dev 153 209.6861^ Test Df^ Deviance 2 level * rotate 151 209.6333 +level:rotate 2 0.
Again, an analysis of deviance indicates that a single treatment odds ratio is sufficient (1-pchisq(.05,2)= 0.98 ). 2e) Summary of the evidence for a learning effect from the continuation ratio logit model: Modeling time 1 responses, the odds of taking longer to complete the task (relative to less time, for any cutpoint) are 1/exp( -2.0075675) = 7.4 times higher for the group that sees the rotated maze, relative to the control group. At time 8, the odds ratio is 1/exp(-0.8071621)= 2.2. The two groups perform more similarly at time 8 than at time 1. This is indirect evidence for a learning effect. We don’t know that either group improved by completing the task in fewer tries – we just know that the gap between the groups is smaller. As in question 1, we can use the EDA to support the idea that the group seeing the rotated map improved in its peformance.
