Lecture 15 (lab 4) —February, January 7, 2006
What was covered?
- R topics
- generic functions in R
- using the apply function to apply functions separately to the rows or columns of matrices
- the surface plotting functions persp (base graphics) and wireframe (lattice graphics)
- the matrix function for creating a matrix from a vector
- the contour graph function contour
- the matrix functions solve and diag
- graphing loglikelihood functions of two variables—surface plots and contour plots
- obtaining maximum likelihood estimates for loglikelihoods of two parameters
- using the Hessian matrix to obtain standard errors for Wald confidence intervals
- generating and interpreting profile likelihood confidence intervals, both joint and marginal
- the Pearson chi-square and G2 goodness of fit tests for models fit to count data
R functions and commands demonstrated
- apply is used to apply a function separately to the rows or columns of a matrix
- contour is used to produce contour plots
- diag extracts the diagonal entries from a matrix
- dim is used to obtain the dimension (number of rows and columns) of a matrix or data frame
- expand.grid is used to create all possible pairings of the elements of two vectors
- matrix is used to convert a vector into a matrix. By default the matrix is filled column by column
- names is used to display or change the names of the elements of a vector
- persp is the 3-dimensional surface plotting function of R
- solve when given a matrix as its argument returns the inverse of the matrix
- trellis.par.get is used to obtain settings for lattice graphics
- trellis.par.set is used to change settings for lattice graphics
- wireframe is the lattice package equivalent to the base graphics function persp. Used to produce 3-dimensional surface plots
- \n is a return code used inside text strings to cause a line break in the printed text
R function options
- colorkey= (argument to wireframe) takes on values TRUE or FALSE. Setting it to FALSE suppresses the display of the color code key generated when the drape argument is set to TRUE.
- drape= (argument to wireframe) takes on values TRUE or FALSE and can be used to add color shading to a surface plot
- nlevels= (argument to contour) is used to specify the number of contour lines to display
- nrow= (argument to matrix) specifies the number of rows for the matrix
- par.settings= (argument to wireframe) used for changing default lattice parameter settings for the current graph only
- phi= (argument to persp) changes the viewing direction. phi specifies the colatitude
- scales= (argument to wireframe) can be used, e.g., to display tick marks on the axes
- screen= (argument to wireframe) can be used to specify a different viewing direction of the surface plot
- theta= (argument to persp) changes the viewing direction. theta specifies the azimuthal direction
- ticktype= (argument to persp) adds tick marks to 3-D surface plots
- type= (argument to plot) controls what is plotted. We used type='l' in this exercise to plot a curve
R packages used
Creating a negative binomial loglikelihood function
- We continue with the data set from last time. Our goal this time is to fit a negative binomial distribution to these data. Since the negative binomial is a two-parameter distribution, this will introduce some new complications.
> #generate data
> num.stems<-c(6,8,9,6,6,2,5,3,1,4)
> aphids<-0:9
> aphid.data<-rep(aphids,num.stems)
- If we paralleled exactly our Poisson loglikelihood function from last time, then the negative binomial loglikelihood function should be the following.
> #negative binomial loglikelihood
> NBfunc1<-function(mu,theta) sum(log(dnbinom(aphid.data,mu=mu,size=theta)))
As you can see the only difference from the Poisson function is that we need to specify two arguments, the mean and the dispersion parameter.
- In truth though it will be more convenient if we use a slightly different version of this function. In this version instead of specifying the two arguments as separate variables, we specify them as components of a vector.
> #alternative vector version
> NBvec.pos<-function(p) sum(log(dnbinom(aphid.data,mu=p[1],size=p[2])))
- The next few lines contrast how the two versions of the negative binomial loglikelihood function differ.
> #first function
> NBfunc1(3,4)
[1] -116.1379
> #incorrect use of 2nd function with two separate arguments
> NBvec.pos(3,4)
Error in NBvec.pos(3, 4) : unused argument(s) ( ...)
> #correct usage in which arguments are entered as a vector
> NBvec.pos(c(3,4))
[1] -116.1379
- For use with nlm recall we need the negative loglikelihood function. I use the function NVvec.pos created above to define this function.
> #negative loglikelihood for nlm
> NBvec.neg<-function(p) -NBvec.pos(p)
Graphing the negative binomial loglikelihood
Two-dimensional profiles
- A complete picture of the negative binomial loglikelihood requires three dimensions. We can though plot two-dimensional profile loglikelihoods by holding one parameter fixed. Since μ is the mean of the negative binomial distribution, we suspect that its maximum likelihood estimate will be the sample mean. I plot profile likelihoods with μ held at the sample mean as well as at nearby values of μ. The following code generates four different profiles. I use the type='l' argument in the plot statement to plot line segments rather than the default points.
> plot(seq(0.1,10,.1),sapply(seq(0.1,10,.1), function(x) NBvec.pos(c(3.46,x))), xlab=expression(theta), type='l', ylab='loglikelihood')
> lines(seq(0.1,10,.1), sapply(seq(0.1,10,.1),function(x) NBvec.pos(c(4.5,x))), col=2)
> lines(seq(0.1,10,.1), sapply(seq(0.1,10,.1),function(x) NBvec.pos(c(3,x))), col=4)
> lines(seq(0.1,10,.1), sapply(seq(0.1,10,.1),function(x) NBvec.pos(c(2.5,x))), col=3)
> legend(6.5,-130, c(expression(mu==3.46), expression(mu==4.5), expression(mu==3), expression(mu==2.5)),
col=c(1,2,4,3),lty=rep(1,4),bty='n')
- In the plot and lines statements I repeatedly make use of what's called a generic function. The syntax I use for the generic function in the plot statement, e.g., is the following.
function(x) NBvec.pos(c(3.46,x)))
What this construction does is to temporarily turn what is a two-parameter function into a one-parameter function. It permits me to specify a value for one of the parameters, μ, while allowing the other parameter, θ, to remain free.
- The alternative to using a generic function here would have been to create a new function in which μ was explicitly specified, but this would need to be done for each of the values of μ to be plotted. The generic function approach is much more efficient and doesn't needlessly create new objects in the R workspace.
Three-dimensional surface plots
- Since the likelihood is a function of two variables, we need three dimensions to effectively display it. Just as we need a sequence of values to plot a curve in two-dimensions, we need a two-dimensional grid of values to plot a surface. R provides a number of functions for this purpose of which expand.grid is arguably the most useful.
- The expand.grid function takes as its arguments two vectors and generates all possible combinations of their components. I illustrate its use in a simple example below.
> #Generating points in the plane
> expand.grid(1:3,8:11)
Var1 Var2
1 1 8
2 2 8
3 3 8
4 1 9
5 2 9
6 3 9
7 1 10
8 2 10
9 3 10
10 1 11
11 2 11
12 3 11
Observe that the elements of the first vector are cycled through first.
- Next I generate the z-coordinate. One way to do this is to just create an ordinary function of two variables and use the columns of the output from expand.grid as input to this function. Here's an example.
> #create the z-coordinate
> test2<-function(x,y) 2*x+y
> expand.grid(1:3,8:11)->yuk
> test2(yuk[,1],yuk[,2])
[1] 10 12 14 11 13 15 12 14 16 13 15 17
- A second way is to use a function that expects a vector as input and then use the apply function to apply the function to the output from expand.grid one row at a time. The next line of code does it this way.
> #another way of doing it with vector argument and apply
> apply(yuk,1,function(x) 2*x[1]+x[2])
1 2 3 4 5 6 7 8 9 10 11 12
10 12 14 11 13 15 12 14 16 13 15 17
- The apply function can be used to apply a function to either the rows, 1, or the columns, 2, of a matrix. The first argument is the matrix to which the function should be applied, the second argument is either 1 or 2 indicating rows or columns, and the third argument is the function written in such a way that it can accept a vector as input and return a single number in return. In the code above I rewrite the test2 function above in vector form treating it as a generic function so that I am able to create it on the fly.
- With that background I create a grid for our loglikelihood function and then evaluate the function on that grid. Based on our profile likelihood plot above I generate values in the intervals 2 ≤ μ ≤ 5 and 1 ≤ θ ≤ 4 in increments of 0.1.
> #create x-y grid for plotting loglikelihood
> g<-expand.grid(mu=seq(2,5,.1),theta=seq(1,4,.1))
> #examine first ten values
> g[1:10,]
mu theta
1 2.0 1
2 2.1 1
3 2.2 1
4 2.3 1
5 2.4 1
6 2.5 1
7 2.6 1
8 2.7 1
9 2.8 1
10 2.9 1
> #create and append the z-coordinate, the loglikelihood
g$z<-apply(g,1,NBvec.pos)
dim(g)
[1] 961 3
The dim function returns the dimensions of our matrix. We see that there are 961 rows and 3 columns.
Method 1: Using the persp function
- The persp function of R can be used to make a serviceable surface graph. A portion of the help screen and an excerpt from the details section are shown below.
> #3-dimensional surface graph. Look at help on persp
>?persp
Fig. 2 Portions of persp help screen
- Observe that persp expects the z-coordinates to be arranged in a matrix whose dimensions match the dimensions of the grid we've created. The rows of the matrix should correspond to the x-values. We can obtain this format using the matrix command. The first argument to matrix is the vector that is to be converted to a matrix. I use the nrow argument to specify how many distinct x-values there are (equal to the number of rows in the matrix). By default the matrix function fills the matrix one column at a time, which is what we want for our data.
> #we construct a matrix of z-values to match the grid
> z.matrix<-matrix(g$z,nrow=length(seq(2,5,.1)))
- I next create the 3-dimensional plot. The persp function does not support mathematical characters so I use the English equivalents for the Greek letters to identify the parameters.
> persp(seq(2,5,.1), seq(1,4,.1), z.matrix, xlab="mu", ylab="theta", zlab='loglikelihood')
- By default tick marks are not displayed. I can add them with the ticktype='detailed' argument.
> #add ticks using ticktype option
> persp(seq(2,5,.1), seq(1,4,.1), z.matrix, xlab="mu", ylab="theta", zlab='loglikelihood', ticktype='detailed'
- The spatial orientation of the graph can be changed by specifying values for the viewing angle arguments theta and phi. The finished result is shown in Fig. 3.
> #change viewing position
> persp(seq(2,5,.1),seq(1,4,.1), z.matrix,xlab="mu", ylab="theta", zlab='loglikelihood', ticktype='detailed', theta=30, phi=30)
Method 2: Using lattice graphics
- R possesses a second system of graphics functions defined in the lattice package. Lattice graphics are R's version of the trellis graphics of Splus. We will find lattice graphs to be especially useful when we fit multilevel models later in the course. Lattice graph functions operate very differently from the base graphics functions we've used thus far.
- The lattice graph equivalent to persp is a function called wireframe. First we need to load the lattice package.
> library(lattice)
- Unlike persp, wireframe expects the z-coordinates to fill a vector, not a matrix. To plot the surface we need to use formula notation in which the z-coordinates are on the left side of the formula followed by the formula symbol, ~, followed by the variables defining the x- and y-coordinates separated by an asterisk. wireframe also has a data argument so I can use just the variable names in the formula.
> #wireframe just takes ordinary vectors
> wireframe(z~mu*theta,data=g)
- The default background is gray. We can change this to "transparent" for the current plot only by including the par.settings argument to wireframe as shown below.
> #change background to transparent for current plot only
> wireframe(z~mu*theta, g, par.settings=list(background=list(col="transparent")))
- Alternatively we can change the background setting for all subsequent graphs using the trellis.par.get and trellis.par.set functions. The trellis.par.get function is used to get the current settings which I copy here to a variable. I then change the settings in this copy and use the copy to change the actual settings with the trellis.par.set function. The first argument to trelllis.par.set is the setting to be changed (in quotes) and the second argument is the name of the object that contains the new settings.
> #get current background settings
> back.color<-trellis.par.get("background")
> back.color
$alpha
[1] 1
$col
[1] "#909090"
> #change color to transparent
> back.color$col <-"transparent"
> #set new background settings
> trellis.par.set("background",back.color)
To see all the lattice graphic settings issue specify trellis.par.get(), i.e., the function with no arguments. Be warned that what's displayed is a very long list.
- It's also possible to change the viewpoint with the screen argument. I illustrate its use next.
> #can also get different viewpoint
> wireframe(z~mu*theta, data=g, screen = list(z = 20, x = -70, y = 3)
- Unlike persp, the wireframe function does support math symbols so we can use Greek letters to label the axes. To get tick marks to display we need to use the scales argument in the manner shown below.
> #add tick marks and nice labels
> wireframe(z~mu*theta, data=g, xlab=expression(mu), ylab=expression(theta), zlab='loglikelihood', scales = list(arrows = FALSE))
- To get the z-label to display without running off the screen I use the return code \n inside the text string to tell R to start a new line at this point. Note: this is why in the read.table command path names were specified with forward slashes / or double backslashes \\. Single backslashes are reserved for specifying special formatting within text strings.
> #split z-axis name
> wireframe(z~mu*theta, data=g, xlab=expression(mu), ylab=expression(theta), zlab="log\nlikelihood", scales = list(arrows = FALSE))
- To obtain color shading of the surface include drape=TRUE. Note: the color key that appears to the right of the graph can be suppressed with colorkey=FALSE.
> #add color with drape argument
> wireframe(z~mu*theta, data=g, xlab=expression(mu), ylab=expression(theta), zlab="log\nlikelihood", scales = list(arrows = FALSE), drape=TRUE)
- Finally I reduce the font size of the z-axis label and z-axis tick marks in order that the entire label is displayed. I do this with the par.settings option again. The results are shown in Fig. 4.
> wireframe(z~mu*theta, data=g, xlab=expression(mu), ylab=expression(theta), zlab="log\nlikelihood", scales = list(arrows = FALSE), drape=TRUE, par.settings=list(par.zlab.text=list(cex=.75), axis.text=list(cex=.75)))
- From the graph we can conclude that there is a greater curvature of the likelihood in the μ direction than in the θ direction. In the θ direction it appears that the likelihood consists of a gently sloping ridge. These differences will be reflected in the confidence intervals we obtain for our point estimates. We would expect much wider confidence intervals for θ than for μ.
Method 3: Contour plots
- While we can learn details about the shape of the loglikelihood from 3-D surface plots, we would be hard-pressed to actually guess the maximum likelihood estimate from such a plot. For that we need a contour plot.
- The base graphics function contour can be used to produce contour plots. The specification of the x-, y-, and z-variables matches what we did for persp. Thus we specify the location for the grid lines in the x- and y-directions followed by a matrix of z-values that matches the grid specification. Unlike persp, contour supports math symbols in the plot. The nlevels= argument can be used to change the number of contours that are displayed.
> contour(seq(2,5,.1),seq(1,4,.1), z.matrix, xlab=expression(mu), ylab=expression(theta))
> #increase the number of contours
> contour(seq(2,5,.1), seq(1,4,.1), z.matrix, xlab=expression(mu), ylab=expression(theta), nlevels=30)
- From the contour plot we can easily read off that the MLE for μ is approximately 3.5 while that for θ is roughly 2.5. Also the difference in curvature in the loglikelihood surface in the μ and θ directions is readily apparent.
Obtaining the MLEs
- I use the nlm function and the negative loglikelihood function created earlier to numerically approximate the MLEs. I use our graphical work to come up with initial guesses for
and 
.
> nlm(NBvec.neg,c(3,4),hessian=TRUE)
$minimum
[1] 114.7009
$estimate
[1] 3.459995 2.645024
$gradient
[1] -2.252383e-05 1.319529e-05
$hessian
[,1] [,2]
[1,] 6.2589599224 0.0002248452
[2,] 0.0002248452 0.9535298199
$code
[1] 1
$iterations
[1] 8
- From the output we see the algorithm converged (code = 1) and both components of the gradient are close to zero. The reported estimate and the minimum value of the negative loglikelihood are close to what we anticipated based on the graphs.
Wald confidence intervals
- To calculate the Wald confidence intervals we need to invert the hessian matrix. Recall that since the hessian reported by R is based on the negative loglikelihood, it is in fact identical to the observed information matrix. The inverse of the information matrix is the variance-covariance matrix of the parameter estimates. The solve function of R when given a matrix returns the inverse of that matrix.
> #save results and hessian
> nlm(NBvec.neg,c(3,4),hessian=TRUE)->out
> #invert the hessian to obtain the covariance matrix
> solve(out$hessian)
[,1] [,2]
[1,] 1.597710e-01 -3.767448e-05
[2,] -3.767448e-05 1.048735e+00
- The diagonal entries are the variances of and respectively. We can extract them from the matrix with the diag function.
> diag(solve(out$hessian))
[1] 0.1597710 1.0487360
- The point estimates are in out$estimate.
> out$estimate
[1] 3.459995 2.645024
- I give the components more informative names using the names function and construct the Wald confidence intervals using the formulas
> names(out$estimate)<-c('mu','theta')
> low<-out$estimate-qnorm(.975)*sqrt(diag(solve(out$hessian)))
> high<-out$estimate+qnorm(.975)*sqrt(diag(solve(out$hessian)))
> cbind(low,high)->wald
> wald
low high
mu 2.6765704 4.243419
theta 0.6378678 4.652180
Profile likelihood confidence intervals
- Recall from last time (Lecture 11) that a (1–α)% profile likelihood confidence interval is defined by the value of the loglikelihood for a test model that would just cause the likelihood ratio statistic to become statistically significant at level α. When there is more than one parameter it turns out there are two distinct kinds of profile likelihood confidence intervals.
- One is the joint profile likelihood confidence interval that is a 95% confidence interval for both μ and θ simultaneously. It is based on the following likelihood ratio statistic.
- The second is a marginal profile confidence interval that is a 95% confidence interval for one parameter while the other is held fixed at its maximum likelihood estimate. This is typically what one means when one says a profile likelihood confidence interval. Here we get intervals for μ and θ separately.
- Observe that the degrees of freedom for the chi-squared statistic are different for the two kinds of confidence intervals.
- We can display these confidence intervals using R's contour function. In order for them to display completely I first recreate the grid specifying a larger range of values for θ.
> #extend graph and create new z-values
> g<-expand.grid(mu=seq(2,5,.1), theta=seq(1,9,.1))
> g$z<-apply(g,1,NBvec.pos)
> z.matrix<-matrix(g$z, nrow=length(seq(2,5,.1)))
> #limits for marginal profile likelihood CI
> lower.limit1<- -out$minimum-.5*qchisq(.95,1)
> #add limit for joint CI
> lower.limit2<- -out$minimum-.5*qchisq(.95,2)
> contour(seq(2,5,.1), seq(1,9,.1), z.matrix, xlab=expression(mu), ylab=expression(theta), levels=c(lower.limit2, lower.limit1), lty=c(3,1), col=c(4,1))
> points(out$estimate[1], out$estimate[2], pch=16, col=2, cex=1.5)
> legend(2, 8.5, c('95% joint CI', '95% marginal CI'), lty=c(3,1), col=c(4,1), bty='n', cex=c(.8,.8))
- The outer contour is the 95% joint confidence interval for both μ and θ. We are 95% confident that the true values for both μ and θ lie within this contour.
- The inner contour defines the outer limits for the marginal profile confidence intervals. To see how this works, we can use the plkhci function from the Bhat package to calculate these intervals and then display them on the graph.
> library(Bhat)
- Recall that the plkhci function requires a peculiar list argument in which the names of the parameters are given as well as estimates for the MLEs and bounds on the confidence intervals.
> x.in<-list(label=c('mu', 'theta'), est=out$estimate, low=c(2.5,1), high=c(4.8,9))
> plkhci(x.in, NBvec.neg,'mu')->profile.mu
> plkhci(x.in, NBvec.neg,'theta')->profile.theta
> rbind(profile.mu, profile.theta)->profile.ci
> colnames(profile.ci)<-c('low', 'high')
> profile.ci
low high
profile.mu 2.748964 4.371168
profile.theta 1.321214 6.465232
- I add these to the graph as horizontal and vertical lines.
> #show boundaries of marginal CI
> abline(v=c(profile.mu),col=3,lty=2)
> abline(h=c(profile.theta),col=2,lty=2)
Goodness of fit
- I produce a graph that compares the fit of the negative binomial and Poisson models (the latter was estimated last time).
> dnbinom(0:9,mu=out$estimate[1],size=out$estimate[2])
[1] 0.10943985 0.16405655 0.16945422 0.14869883 0.11893284 0.08958118
[7] 0.06468935 0.045278220.03093792 0.02073880
> dnbinom(0:9,mu=out$estimate[1],size=out$estimate[2])*50
[1] 5.471993 8.202827 8.472711 7.434941 5.946642 4.479059 3.234468 2.263911
[9] 1.546896 1.036940
> #anything left over?
> sum(dnbinom(0:9,mu=out$estimate[1], size=out$estimate[2])*50)
[1] 48.09039
- I lump counts 9 and above in the last category
> expprobs<-c(dnbinom(0:8,mu=out$estimate[1], size=out$estimate[2]), 1-pnbinom(8, mu=out$estimate[1], size=out$estimate[2]))
> expprobs*50->Ei
> table(aphid.data)->Oi
- I rename the last category so that it indicates 9 and above.
> names(Oi)[length(Oi)]<-"9+"
- I next produce a bar graph that compares the fit of the Poisson and negative binomial models.
> coords<-barplot(Oi, ylim=c(0,11))
> points(coords,Ei,col=2, pch=16, cex=1.2)
> lines(coords,Ei,col=2)
> #add poisson
> poisprobs<-c(dpois(0:8, lambda=3.46), 1-ppois(8, 3.46))
> poisprobs*50->ei.pois
> ei.pois
> points(coords, ei.pois, col=4, pch=15, cex=1.2)
> lines(coords, ei.pois, col=4)
> legend(coords[7],11, c('Poisson', 'negative binomial'), col=c(4,2), lty=c(1,1), pch=c(15,16), bty='n', cex=c(.9,.9))
> box()
- Graphically the fit of the negative binomial looks far better than that of the Poisson. I next test this formally. If we examine the expected counts we see that there are a few below 5.
> Ei
[1] 5.471993 8.202827 8.472711 7.434941 5.946642 4.479059 3.234468
[8] 2.263911 1.546896 2.946552
- Recall that if the expected counts are small (5 being a guidline for smallness) that the chi-squared distribution of the Pearson statistic becomes questionable. To deal with this I combine categories. On way to do this would be to leave the first five categories alone, combine the next two, and combine the last three.
> Ei.merge<-c(Ei[1:5],sum(Ei[6:7]),sum(Ei[8:10]))
> Ei.merge
[1] 5.471993 8.202827 8.472711 7.434941 5.946642 7.713526 6.757360
- I repeat this with the observed counts and carry out the test. The degrees of freedom for the chi-squared statistic is m – 1 – p where m is the number of categories (after merging) and p is the number of estimated parameters (two here corresponding to μ and θ).
> Oi.merge<-c(Oi[1:5],sum(Oi[6:7]),sum(Oi[8:10]))
> Pearson<-sum((Oi.merge-Ei.merge)^2/Ei.merge)
> Pearson
[1] 0.6607193
> Pearson.p<-1-pchisq(Pearson,length(Oi.merge)-1-2)
> Pearson.p
[1] 0.9560832
- The result is not even close to being significant. I next carry out the G2 test and compare the results.
> #G2 test
> G2<-2*sum(Oi.merge*log(Oi.merge/Ei.merge))
> G2
[1] 0.6675889
- The results match that of the Pearson test. What is the critical value for this test?
> qchisq(.95,length(Oi.merge)-1-2)
[1] 9.487729
- So we should reject the null hypothesis and conclude the lack of fit is significant if the observed value of our test statistic, Χ2 or G2, is greater than 9.49. That is clearly not the case here. Observe that the critical value is roughly twice the degrees of freedom.
> length(Oi.merge)-1-2
[1] 4
- What would happen if we didn't merge cells?
> Pearson<-sum((Oi-Ei)^2/Ei)
> Pearson.p<-1-pchisq(Pearson,length(Oi)-1-2)
> Pearson
[1] 3.511331
> Pearson.p
[1] 0.8340238
- Of course the chi-squared distribution of the test statistic is suspect here. We can obtain a better estimate of the p-value by using a randomization test.
> chisq.test(Oi,p=expprobs,simulate.p.value=TRUE)
Chi-squared test for given probabilities with simulated p-value
(based on 2000 replicates)
data: Oi
X-squared = 3.5113, df = NA, p-value = 0.9455
- The result is nearly identical to what we obtained with the merged categories.
- A skeptic might argue that by merging the expected probabilities into the last category we unjustly favored our model (because it makes the expected results more closely resemble the observed data). To satisfy such a skeptic I create a new category X = 10+ for both the observed and expected counts and redo the randomization test. The frequency of this new category for the observed data is 0.
> new.oi<-c(Oi,0)
> new.expprob<-c(dnbinom(0:9,mu=out$estimate[1],
size=out$estimate[2]),1-pnbinom(9,mu=out$estimate[1],
size=out$estimate[2]))
> new.expprob
[1] 0.10943985 0.16405655 0.16945422 0.14869883 0.11893284 0.08958118
[7] 0.06468935 0.04527822 0.03093792 0.02073880 0.03819224
> new.expprob*50
[1] 5.471993 8.202827 8.472711 7.434941 5.946642 4.479059 3.234468
[8] 2.263911 1.546896 1.036940 1.909612
> chisq.test(new.oi,p=new.expprob,simulate.p.value=TRUE,B=2000)
Chi-squared test for given probabilities with simulated p-value
(based on 2000 replicates)
data: new.oi
X-squared = 13.5113, df = NA, p-value = 0.1844
- So even after creating a separate category in which the observed count is zero we still fail to find a significant lack of fit.
Cited References
- Krebs, Charles. 1999. Ecological Methodology. Addison-Wesley: Menlo Park, CA.
Course Home Page
