title: Plotting a gaussian mixture with pdf values >> 1 in MATLAB

text: From a given dataset X, I learn a 7-component gaussian mixture model using matlab's gmdistribution.fit:               model = gmdistribution.fit(X,7,'Start', gmm_init,'Options',options, 'Regularize', 1e-05)      However, if I want to create a contour-plot using matlab's ezcontour, the result is that a lot of areas are actually blanked out in the plot (matlab's plot data cursor tells me these entries are "NaN"): ![enter image description here](http://i.stack.imgur.com/nFXHr.png) If I look at the pdf of my model evaluated at X's data points, I get that a lot of pdf values are >> 1:                >>PDF = pdf(model, X);      >>max(PDF)           ans =              2.806149789376303e+13           >> mean(PDF)           ans =              6.587701947551628e+12      From other questions I have learned that in general, a continuous PDF may indeed be larger than 1. However, since in my case the average pdf value is some orders of magnitude higher than the smallest one, it seems that ezcontour has some difficulties plotting it. My dataset X contains some very dense regions. Also, some mixture components only contain very few components (e.g. ~1% of all data points). In a different post, it is suggested that I look at the (co)variance of my components to get an idea about the density for a given component's assigned data points. My questions are:   1. Is the model I get from gmdistribution.fit of an inferior quality because the PDF has some really crazy high values in the dense,"middle" regions of my dataset?   2. What can I do to plot this gaussian mixture model anyway (e.g. with something else than ezcontour?)


title: Leave-one-out cross-validation in LiblineaR

text: So I've been encountering this weird issue with LiblineaR where I get different results from just running LiblineaR with cross set to the size of my dataset and actually looping through the dataset, choosing each entry for testing and training on the rest, even though it seems like these should be exactly the same. Has anyone else dealt with this? (This is LiblineaR's LOO cross-validation)               LiblineaR(data=data,labels=yTrain,cost=co,cross=nrow(data))      (This is the code I've been using, so that I can see which specific works were misclassified and save the weights of the training model if necessary)               for (i in 1:nrow(data)){     x = data[,2:dataCols]     y = factor(data[,1])          xTrain = x[-i,]     xTest = x[i,]     yTrain = y[-i]     yTest = y[i]          m=LiblineaR(data=xTrain, labels=yTrain, cost=co)     p=predict(m, xTest)}


title: Estimating correlated parameters with multi-level model

text: I would like to estimate a multi level model in Stata or R (using lmer) where the first level coefficients are the same for all observations, but the coefficients within observation are correlated. An example would look something like this: $$Y_i = \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + ... + \varepsilon_{0i}$$ $$\beta_1=\gamma_1 z_1 + \gamma_2 z_2 + \varepsilon_{1}$$ $$\beta_2=\gamma_1 z_1 + \gamma_3 z_3 + \varepsilon_{2}$$ $$\beta_3=\gamma_2 z_2 + \gamma_3 z_3 + \varepsilon_{3}$$ and so on, with equations for each beta. Clearly, I'd make a distributional assumption for the $\varepsilon$'s... like $\varepsilon \sim N(0,\sigma^2)$ The x variables vary by observation, but the z variables do not vary between observations. Thus, the parameters $\gamma$ and $\beta$ are also the same for all observations. This differs from most hierarchical models I have seen in that parameters are related within an observation, rather than depending on observation-level characteristics. As a specific application, consider a model where the dependent variable $Y$ is a student's test scores. The x variables are measures of performance in previous classes, and the $z$ variables are characteristics of those classes. Students have taken the same set of classes, but there are few students in each class, so I'd like to pool estimation of the coefficients $\beta$. Because the classes have similar characteristics, there may be far fewer $\gamma$ parameters than $\beta$ parameters, and pooling estimates to those lower level class characteristics may yield more precise estimates of $\beta$ than estimation without the 2nd level model. At the same time, I'd like estimates of the $\beta$ parameters, so substituting in and estimating y as a function of $\gamma$ and x only gets me half way there. What is the best way to estimate this type of model? I typically program in R, Stata and Python.


title: Distance to a regression line, and degrees of freedom

text: How do you estimate degrees of freedoms for derived measurements? I want to assess the significance of the distance of an independent data point to a regression line. I can easily calculate the (vertical) distance between the data point and the regression line, and I get the uncertainty of the distance from the uncertainties of slope and intercept of the linear regression via Gaussian error propagation. However, what are the degrees of freedom? The linear regression line has been calculated from `n` data points, thus its degrees of freedom is `n-2`. The additional measurement is independent, so I get another degree of freedom, bringing the total to `n-1`? Also, should I estimate the uncertainty of the independent measurement using the variance of the residuals of the regression, since the measurement process is the same for both the data that went into the fit and the independent data point? I guess this would reduce the degrees of freedom again?


title: Disadvantages of negbin regression

text: Does anybody know if there're common known disadvantages of a negbin regression? In my opinion it seems to fit every problem pretty good (measured with the estimated dispersionparameter). So why not always use it?


title: Fastest way to compare ROC curves

text: I have a set of true positive (TP) values which are used to train a model. I am using 5-fold cross validation to train my model (i.e. split my true positives into 5, use 4/5ths for training and 1/5th for testing) I repeat this using different 1/5ths as the test set. For each run, I have a large set of mixed true positives / true negatives which I use my trained model to attempt to classify. I then obtain an ROC curve. This is done for each run of the cross validation (i.e. I end up with 5 ROC curves) I then average the AUC and return it. ### My problem: I have two methods of classification: call them method A and method B. for each method, I get 5 ROC curves. How can I determine which method gives me a better ROC if I have more than one ROC for each? I know computing the AUC and averaging for each method, then comparing averaged AUCs is NOT a good approach. Note: I have more than 1 model (roughly 120). I just explained in terms of one model for simplicity. So I have 120 models, each one having classified the data using method A and method B, and for each method A & B there are 5 ROCs from cross validation. ## Edit * * * My problem more specifically is that I have >100 sets of sequences, and for each set I construct a position weight matrix, which I then use to score against all sets merged together. I have several scoring schemes so I'd like to determine which ones give me the best classification. For this, I use cross validation: split my data into 5 for each set, train my pwm with 4/5ths of the data and test it on 1/5th. Pool the results from 5 runs, and plot an AUC.


title: Estimating the confidence interval for each half of a population

text: We are looking at how a change in the environment (confinement level) alters which of two behaviors our model organism chooses. We perform several trials at each confinement level. Each of these trials elicits one of the two possible behaviors. From this we get a fraction of trials which exhibited each behavior under each of the environmental conditions. I want to calculate the confidence interval for each of the fraction estimates under each environmental condition. So an example would be at a confinement level of 10 (on some arbitrary scale) behavior one is used in 23 of 31 trials, and behavior two is used in the other 8. My guess is that the proper thing to do it to use r's 1−sample proportions test with alternate k values for each pair of fractions (as in: http://www.r-tutor.com/elementary-statistics/interval-estimation/interval- estimate-population-proportion). So for our example, prop.test(23, 31) and prop.test(8, 31) giving confidence intervals of (0.55,0.87) and (0.12, 0.44). Or is there some reason not to use a 1-sample proportion test on each member of the fraction pair? If the answer is just, "nope, that's fine", well... great!


title: Which is largest, of a bunch of normally distributed random variables?

text: I have random variables $X_0,X_1,\dots,X_n$. $X_0$ has a normal distribution with mean $\mu>0$ and variance $1$. The $X_1,\dots,X_n$ rvs are normally distributed with mean $0$ and variance $1$. Everything is mutually independent. Let $E$ denote the event that $X_0$ is the largest of these, i.e., $X_0 > \max(X_1,\dots,X_n)$. I want to calculate or estimate $\Pr[E]$. I'm looking for an expression for $\Pr[E]$, as a function of $\mu,n$, or a reasonable estimate or approximation for $\Pr[E]$. In my application, $n$ is fixed ($n=61$) and I want to find the smallest value for $\mu$ that makes $\Pr[E] \ge 0.99$, but I'm curious about the general question as well.


title: Why are there no one-inflated count data models?

text: I am working on zero-inflated count data models using the `pscl` package. I am just wondering why there is no development of models for one-inflated count data models! Also why there is no development of bimodal, say zero- and-2-inflated, count data models! Once I generated one-inflated Poisson data and found that neither the `glm` with `family=poisson` model nor the negative binomial (`glm.nb`) model was good enough to fit the data well. If any one can shed some light on my thought, eccentric though it might be, it would be very helpful for me.


title: How can I generate correlated timeseries made up of 0s and 1s?

text: I want to generate series of 0s and 1s that exhibit some clustering. By this I mean that 1s and 0s should occur together. So I envisage series of 0s and 1s that will exhibit similar clustering of these elements, and not just random series of 0s and 1s. In essence, for a single time series, I would go about that by thresholding a Markov chain with a 2x2 transition matrix, with some stochastic added on to it. Now, I'm not too certain on how to do this, but since I would like to produce several of these series, I was wondering whether there's something straightforward that I have missed. I plan to use these series to simulate data availability (0 or 1) in a data acquisition system and to do some Monte Carlo simulations of how this affects what we can do with the data. In order for the above simulations to be realistic, I would like to fit real observations to this model, so that the simulated time series share temporal correlation with the data. I would initially do this by calculating lag autocorrelations of both series and tweaking model parameters until I get something that resembles my observations, but unsure whether this is the best way. Thanks!


title: Cross validation training support vector machine

text: Can somebody explain me how to perform cross validation while training Support vector machine via some example?


title: Quantifying effect of a categorical variable in time series analysis

text: I have a dependent variable ($y$ axis in the pictured graph) trending over time ($x$ axis). I also have a categorical variable, which was constant until one point in time, then changed to another value, with which value it has remained constant since; that point in time is indicated in the picture by a break in the pictured graph. What test can I use to find the effect size of the categorical variable? ![trend](http://i.stack.imgur.com/Xf88J.png)


title: How to apply an AR(MA) model to a prewhitened signal?

text: I have two (vehicle velocity) signals that should consist of similar "latent" drivers, but have different autocorrelation structures. The driver-signals are quite nasty statistically, so I'm not attempting to model them. I can get quite nice results by prewhitening the signals using AR(1)-residuals, but these are very difficult to interpret in "real world terms" (ie. velocities). So what I'd like to do is to prewhiten one of the signals and then add the AR-model of the other signal to this, so that I'd have two signals with same autocorrelation structures. It may be that there is a very simple method for doing this, but unfortunately I haven't found one, or it maybe impossible. I guess it should be sort of an inverse of the Yule-Walker method. One also that is quite close is to use arima.sim with innovations, but with the difference that I don't have innovations, but residuals.


title: Estimating covariate effects on Likert-type response variable

text: I am working on a problem where we elicit medical residents perceptions regarding the importance of specific core medical competencies (professionalism, collaboration, medical expertise, etc.) on a 5 point Likert scale. For example, "In your current residency role how important is professionalism in your day-to-day activities". We elicited responses from two schools (School1, School2) and from residents in 5 types of post-graduate residency programs (Family Medicine, Psych, Surgery, Internal, Other). We have data on approximately 800 respondents. Our hypothesis is that the culture surrounding the specific schools in the sample, and the specific residency programs under question may influence the self-reported importance of these roles. Obviously, the response data live on 5 points. And in our case, the responses are fairly skewed to the right, with many residents finding these core roles important. I have graphically investigated the relationship between Perception (5 point Likert Scale) and my two independent variables (School, Residency Program). I have used mean plots, also considered mosaic plots. I have a sense of what is going on qualitatively but want estimated covariate effects, p-values, etc...corresponding to the impact of School, Residency and their interaction on self-reported importance. Ideally, I could run this as a 2-way ANOVA (or linear regression model). However, diagnostic investigation of residuals suggests an issue of non- normality. I considered the proportional odds logistic regression model next; however, in many instances the assumption of "proportional odds" was not satisfied. From there I moved on to the multinomial logistic regression model; however, given the distribution of the response data (skewed to the right), very little information exists in cells {1}, {2}, so estimated covariate effects (relative to these response levels) appear to be estimated without a great deal of precision. Next, I considered non-parametric tests...extensions of Kruskal Wallis test from one-way layouts to two-way (and higher order) layouts. I found the following: http://www.tandfonline.com/doi/abs/10.1207/s15327906mbr1503_4#.Ud1uCdpzY3E However, I could not find a statistical package to implement such a method. Does anyone know of one? I wouldn't expect these tests to elucidate the size of the covariate effect; however, they may at least be able to derive a valid overall p-value (something like the global F-test in regression/ANOVA)...or perhaps even more valuable LRT style p-values corresponding to each covariate effect (i.e. School, Program, School*Program). I imagine this is a fairly common problem in applied research settings...just wonder what other people have done in the past? Thanks Chris


title: Change points in the mean and in the variance using R

text: Hi I am trying to check the existence of change points in the mean and in the variance using R. My data is not normally distributed, so I want to use the cusum test. Also I want the significance level to be 0.05. I am trying to use the package "changepoints", but for change points in the mean I couldn't find a way of choosing the significance level because it says the asymptotic penalty is not available with cusum! For the variance it lets me choose the significance level using: `single.var.css(x,penalty="Asymptotic",value=0.05,class=TRUE,param.estimates=TRUE)` however it always detects a change point. Sometimes it looks like a real changepoint, but when the series doesn't seem to have one it just detects a change point in the end of the series. Anyone knows a way of calculating change points in R using the cusum test with a significance level of 0.05? Thanks for helping!


title: The search for statistical dispersion before tests… and multiple testing issues

text: After finishing a data collection season, there's a huge dataset containing n>>1 variables (columns) and a few different observations divided by group (rows; let's say cases and controls). The task is to find variables with "LARGE" differences between cases and controls, but before starting any test, there's a special feature to consider about the dataset: Descriptive statistics show that most of the variables have very small dispersion (i.e., all observations are very similar, so standard dev. ≈ 0)..., although a few of them (approx. 1% of all variables) do show higher SD values. Then, although differences between cases and controls could -in principle- be detected across many variables, the fact that the majority of these variables do not show large statistical dispersion makes them unatractive. That is to say, small between-group differences would not be interesting, although they were statistically significant... So: Should only variables with high dispersion be chosen to compare groups in statistical tests? That is to say: is it tricky, if one chooses only the subset of the top-10 variables with large statistical dispersion values, and then performes group comparisons (in order to avoid doing tests whose results would have no value)? How would Bonferroni work in this case? Would the multiple testing correction take into account only the number of tests performed, with the top-10 subset variables? Any suggestions and/or references on this topic would be appreciated. Thanks in advance.


title: Different ways of making a set of numbers (all between $0$ and $1$) to sum up to $1$

text: I have a set of numbers $S$, and for each $s_i\in S$, $0\lt s_i \lt 1$. I would like to transform them so that they sum up to $1$. An obvious way to do it is to calculate $t_i=\frac{s_i}{\sum_i{s_i}}$, so that $\sum_i{t_i}=1$. I wonder if there is any other way to achieve the same effect, and what are the differences between them.


title: WinBUGS truncated normal distribution

text: I am estimating a stochastic frontier with a mixed model. So far the half normal distribution worked good but I need a truncated normal distribution. It does not work, and I receive the error „Expected collection operator c“. I am using R2WinBUGS and as you can see in the model I have tried OpenBUGS and WinBUGS. Any suggestion? My model looks like:               for (i in 1:N)      {              mu[i] <- alpha + x[i,1]*beta[1] + x[i,2]*beta[2] + x[i,3]*beta[3] +               x[i,4]*beta[4] + x[i,5]*beta[5] + u0[county[i]] + u1[county[i]]*x[i,1] +               u2[county[i]]*x[i,2] + u3[county[i]]*x[i,3] + u4[county[i]]*x[i,4] +               u5[county[i]]*x[i,5] - z[ID[i]]              y[i] ~ dnorm(mu[i], tau)          }               for (i in 1:220) {         z[i] ~ djl.dnorm.trunc(rho,lambda,0,1000)              z[i] ~ dnorm(rho,lambda)T(0,100)  #For openbugs              eff[i] <- exp(-z[i])     }      prior for `rho~dnorm(0,0.027)` I would appreciate your help! Regards, Daniel


title: Interactions between covariates

text: Suppose I have a sample with two covariates whose levels are associated with some response; e.g., the mean height of a plant at 6 months is modeled by two factors: cultivar of plant (say, $A$ versus $B$) and size of pot (e.g., "small" versus "large"). Furthermore, suppose that the truth is that cultivar $A$ in general grows taller than cultivar $B$ irrespective of the size of the pot in which it is planted. However, the pot size has no effect on height at 6 months. It so happens, though, that through random chance or some underlying hidden relationship, twice as many cultivar $A$ plants got the large pot as cultivar $B$, whereas in the small pot, the split is 50/50. Now, if we naively do a $t$-test for $H_0 : \mu_L = \mu_S$ versus $H_a : \mu_L \ne \mu_S$, we would find significance not because the pot size has an effect, but because it happened that more $A$ cultivars were potted in the large pots. Clearly such a test is not appropriate. However, the question I'm faced with is this: is it correct to say that cultivar and pot size have an interaction? Note that pot size doesn't influence the outcome. I know I should be using a 2x2 factorial design, but I have some difficulty interpreting the meaning of the tests: If I compute a confidence interval based on $(\bar x_{AL} - \bar x_{AS}) - (\bar x_{BL} - \bar x_{BS})$, what is the related hypothesis, and is this test accounting for the fact that $A$ gets planted into large pots more frequently than $B$?


title: Clustering on structural variables?

text: I'm working with land surface models. These models basically take a bunch of meteorological forcing data (downward radiation, wind, rain, humidity, etc), and run it through some biogeochemical/physical equations, and spit out a bunch of output data ($CO_2$ update/release, net radiation, evaporation, runoff, etc.). There are some problems with the models, and I'm trying to use clustering methods to split the input variables up into subdomains, so I can look at the responses in the output in a non-linear manner. My question is whether it is valid and/or sensible to cluster on _structural variables_ (e.g. time and space related variables). The input variables, like downward radiation, are _highly_ dependent on these variables (e.g. the amount of sunshine is largely determined by latitude, time of year, and time of day). So I guess it doesn't make much sense to cluster on both the structural and the input variables. But clustering on the structural variables may offer things that clustering on the input variables doesn't - the ability to categorise model output by time of day, for example. Is there any reason why using structural variables should be avoided? Or is it entirely dependent on the model in question?


title: How to determine exact point of tangency?

text: I have fitted my stress-strain data with $y=ax^3+bx^2+cx+d$ and also added tangent lines as shown in figure below. I am interested to see the deviation of the fit from the linearity. I am not aware of any theory that can be used to find the tangent point (the point of deviation of the fit from the linearity). Hence, based on my visual inspection I somewhat decided my tangent point at added vertical lines. This method is of course raising question of reliability. So I am looking for suggestions how I can quantitively determine the points of deviation for all four data sets. Your advice/suggestion would be highly appreciated! ![enter image description here](http://i.stack.imgur.com/2T6xE.png) Question update: The figure shows that there is a linear relationship at the beginning, but at some point (point of deviation from linearity), the linear relationship disappears and non-linear relationship prevails. I'm looking for a way to determine this point of deviation from linearity quantitatively, not as done by visual inspection. How do I do that? (Sorry, not sure what should be the appropriate tag for this post!)


title: Mediation on an interaction with dichotmous mediator

text: We are using Tingely et al's mediation R package for mediation analysis. Our study is a 2x2 design, say A (0 and 1) and B (0 and 1), with dichotomous mediator M and continuous dependent measure Y. Our model is Y = a + b1* A + b2* B + b3* A*B + cX = u where a, b's,and c are coefficient, u error term. We found that (1) the effect of A alone, i.e. the coefficient b1 was significant, (2) the effect of B alone, the coefficient b2 was significant, and (3) the interaction effect of A and B, the coefficient b2 was significant. We wanted to examine how the mediator M mediates the effect in (1), (2) and (3), and were able to conduct the mediation analysis for (1) and (2) by using the moderated mediation feature of the _mediation_ package. However, we ran into a problem when we tried to do the mediation in (3). Would anyone on the forum advise if it is possible, and how, to do the mediation of an interaction effect in an experiment with a 2x2 design?


title: How to distribute a prize among a group of people given their scores?

text: I asked the same question on math.stackexchange and had some responses, but I would also like to hear some input specifically from statisticians and data analysts and I feel that you guys may have something to offer. The question is about finding reasonable ways of dividing a prize among $n$ people in the following situation: To make the example specific, we have $6$ people in total who are going to share a prize of $100$ dollars, and let us denote the amount received by each person $i$ as $q_i$. In addition, each person $i$ is given a score $s_i$, and we can think of $s_i$ as a way of measuring how well person $i$ deserves some portion of the prize. The intuition here is that we would like a higher- scoring person to receive a larger portion of the prize than a lower-scoring person, that is, $q_i\geqslant q_j$ if and only if $s_i\geqslant s_j$. Further, the scores are bounded, so $s_{min} \le s_i \le s_{max}$. The problematic thing here is that $s_i$ can be either negative or positive. For example, $s_1=1.3, s_2=2.1, s_3=-0.8, s_4=-3.7, s_5=0.7, s_6=5.2$. So, what would be the proper ways of dividing the prize given these scores? One interesting answer by @opt suggests to use the so-called Softmax function in the context of neural networks, and it is basically ${\displaystyle p_i=\frac{\exp(s_i)}{\sum^n_j\exp(s_j)}}$, and $\sum^n_ip_i=1$. In other words, $p_i$ would be the portion of the prize that $i$ should receive given her score. I would like to hear your thoughts/opinions on this method. Many thanks.


title: Expectation of binomial distribution

text: From my understanding for a random variable $X \sim Binom (n,p)$, $E(X)=n \times p$, correct? I was looking at an example problem: Every pair of the $n$ people is friends with probability $p$. What is $E(X)$? $X$= number of friendships. $X\sim Binom (\binom{n}{2}, p)$ and $E(X)=\binom{n}{2}p$. This makes perfect sense. However I got confused by the next part: What is $E(Y)$, where $Y$= number of friendship triangles? And the answer is $E(Y)=\binom{n}{3}p^{3}$. However, I don't understand why the $p$ is cubed?


title: Are misses in my data distributed completely at random?

text: I have a large data set containing children's scores on arithmetic tasks with a lot of missing values possibly due to the age of the children? My hypothesis is that the missingness is not completely at random but at random. I need to do a CFA and want to know whether I need to condition on age or not? I know SPSS 17 has a missing values analysis but my data is in R and is quite large. I would like to do a missing values analysis in R but have not been able to find a package that can do this. Does anyone know?


title: Statistics for multi-test replicated correlation analysis

text: I'm analyzing pairwise correlations of time series between two different types of microarrays done for several samples as biological replicates. So, I have M1 number of variables on type 1 array, M2 number of variable on type 2 array, N samples and T time points. For each sample, I calculate M1 x M2 correlation coefficients (Pearson or Spearman) and p-values using T points. Due to batch effect between samples I cannot average measurements between different samples for each time point. My question is which statistics to use to find which pairs of M1 and M2 variables have statistically significant correlation consistent between samples? How to apply multi-test correction? Please help.


title: How to remove all the non parent-child links in a Directed Acyclic Graph?

text: I have some Directed Acyclic Graph trees representing biological features. In my database, I can retrieve all the links between every child and its ancestors: I have all the edges and the nodes. Unfortunately, I have not only the direct links between every child and its parent, but I have also the indirect links between the child and the ancestors. I want to remove all these links between children and ancestors. You can see an example in the following picture: ![enter image description here](http://i.stack.imgur.com/G7o7O.png) From my link list, I want to remove LS1S4 and LS1S3 and LS1S5 (the red ones), and mantain LS1S2, LS1S3, LS1S4 (the black ones). I need an algorithm to do this, I've tried to design many of them but I always failed something... Here are the input data:   * nodes: S1, S2, S3, S4, S5.   * links: LS1S4, LS1S3, LS1S5, LS1S2, LS1S3, LS1S4. How to understand which are the red links and which the black links? Thanks


title: Correctly specifying temporal pseudo-replication in a mixed effects model in R/nlme

text: I would like to correctly specify that my data is temporally pseudo replicated in a mixed effects model in `R` using `lme()` of the `nlme` package. I have done this following the example in Crawley's `R` text book, but I need some confirmation that the fixed and random effects are specified correctly. I have checked previous postings on similar questions and I see conflicting information. Hence, I decided to ask specifically on this forum. Also, included in the model is a covariate controlling for spatial structure, and I would like to confirm that it should be a fixed effect. The structure of my data: The response variable is turnover in species composition (`Tsc`) between two censuses in 26 transects in a tropical forest. Transects are grouped based on their logging status into four levels (unlogged, lightly logged, moderately logged, and heavily logged). There are three time intervals for which turnover was computed (0, 7, 14). Transects of the same logging status occur in the same location, so Latitude is included as a covariate to control for logging status. My interest is in the effect of logging status and time on turnover. The model I constructed following Crawley's R text is as follows:               m1<-lme(Tsc~Logging.status*Time.interval+Latitude,              random=~Time.interval|Transect, data=mydata)      My questions:   * Is it right for Time.interval to be in both the fixed and random effects? How do I correctly specify that each transect has repeated measures for three time intervals?    * Should Latitude be a random effect since I am not really interested in it per say, but only need it to control for spatial structure?


title: Benefits of using QQ-plots over histograms

text: In this comment, Nick Cox wrote: > Binning into classes is an ancient method. While histograms can be useful, > modern statistical software makes it easy as well as advisable to fit > distributions to the raw data. Binning just throws away detail that is > crucial in determining which distributions are plausible. The context of this comment suggests using QQ-plots as an alternative means to evaluate the fit. The statement sounds very plausible, but **I'd like to know about a reliable reference supporting this statement.** Is there some paper which does a more thorough investigation of this fact, beyond a simple “well, this sounds obvious”? Any actual systematic comparisons of results or the likes? I'd also like to see how far this benefit of QQ-plots over histograms can be stretched, to applications other than model fitting. Answers on this question agree that “a QQ-plot […] just tells you that "something is wrong"”. I am thinking about using them as a tool to identify structure in observed data as compared to a null model and wonder whether there exist any established procedures to use QQ-plots (or their underlying data) to not only detect but also describe non-random structure in the observed data. References which include this direction would therefore be particularly useful.


title: Remove outliers and calculate average

text: I am trying to understand (by reverse engineer) the method a program is using to calculate an average value over a population. The program displays a (very very large) image consisting of tiles. Each tile has some overlap with it's neighbors, expressed in pixels. In total there are 8625 tiles, hence 8625 overlap values (considering the x-axis overlap only). The program calculates a value that looks like an average of all the overlap values, but in some cases it is not equal to the simple average. Here are some examples: If all the tiles have a x-overlap of **50 pixels** , the average is of course **50** and the program also reports **50** as the average value. If **4312** tiles have an overlap of **1 pixel** and the rest **4313** have an overlap of **100** , the average is **50.50573913** , but the program says **50.511583**. If **4312** tiles have an overlap of **10** pixels and the rest **4313** have an overlap of **100** , the average is **55.005217391** , but the program says **55.01053**. If **4312** tiles have an overlap of **50** pixels and the rest **4313** have an overlap of **100** , the average is **75.00289855** , but the program says **75.00585**. If **4312** tiles have an overlap of **99 pixels** and the rest **4313** have an overlap of **100** , the average is **99,500058** and the **program agrees** (it also reports 99.500058 as the average). The program outputs a summary file which states the following: **outlier limit='4'** I don't know what that number means or what units it is expressed in though. Assuming that the difference in this results is not due to rounding errors (I am pretty sure it is not), what approach could I follow in order to understand it's logic of removing outliers, or make them weigh less somehow in the calculation of the average?


title: Regression based for example on days of week

text: I need a little help to move in the right direction. It's a long time since I studied any stats and the jargon seems to have changed. Imagine that I have a set of car-related data such as   * Journey time from town A to town B   * Distance from town A to town B   * Engine size   * Driver's shoe size   * Make and model of car   * Day of week I want to predict journey time. I imagine there's a strong correlation between time and distance and probably a weaker one to engine size (and none to shoe size). Presumably multiple regression analysis / ANOVA is the tool to use. But how do I include day of week, since just coding it as Sunday=1, Monday=2 etc feels very wrong? Having used Excel's regression tool, for example, how do I interpret the results? Presumably if R is close to 1 this is good (although if there are many data items it seems as though it can be small yet still be significant). But some sources refer to r-squared which seems to be the SD, so a value close to zero is good. It also shows the t Stat, P-value, F and Significance F, whatever they may be. Can anyone recommend a good reference source?


title: How to present several ANCOVA findings

text: I have six scatterplots (such as the one below) that show an insignificant interaction of the factor (X2) (e.g. gender, age, education etc.) with the predictor (X1) on Y. (X1 is significant.) ![enter image description here](http://i.stack.imgur.com/cq0GH.jpg) The simplest interpretation is that as X1 increases, Y increases for both levels of the factor (in this case, males and females). In other words, gender does not affect the relationship between X1 and Y and the difference in slope is due to sampling error. These two statements apply to all other factors (e.g. age, which has three categories; education, which has three levels etc.). One option is to present the six scatterplots and repeat the above statements as explanation. However, this seems tedious. **Is there a better (and succint way) of presenting these findings?** (I guess this is not a statistical question but it certainly relates to data analysis and interpretation!)


title: Adding noise to a matrix/vector

text: How to add noise to an object, regardless of the amount of dimensions? I would expect to try to use some sort of `apply` to add noise of the type: `runif(1, -0.00001, 0.0001)` for each cell. What I tried before was to have a matrix of r * c, and then creating a matrix with noise, and add them together. It works, but I believe there is a more general solution (and I am interesting in how `apply` could solve this) As a side question, should we add uniform errors or normally distributed errors? EDIT: The reason is that I am trying to do bootstraps, and my distributions are heavily skewed. An overrepresentation of 1's. The problem is that, when I am trying to do cross-validations, sometimes I end up doing an `lm(Y ~ X1, X2)` In which there is no variance (all 1's), and therefore it can't compute and gives an `NA`. However, this is not good for bootstrapping anymore, because these errors are messing up what is "really" going on. I figured adding noise allows the linear model to run, since there will not be 0 variance any more?


title: including interactions in a mixed linear model in R (lme)

text: I'm trying to test for the effect of soil moisture on transpiration rates. I have plot-level data for 18 plots in 6 different stands of trees (3 plots x 6 stands). I want to treat "stand" as a random effect. There are reasons (genetic similarities within a stand, differences in elevation, species composition, etc.) to think that stands may differ in transpiration rates. Moisture varies by stand as well. The simple model in made in R is: `model<-lme(transpiration~moisture,random=~1|stand,data=data)` Does this satisfy it? Or should I do: `model<-lme(transpiration~moisture+moisture*stand,random=~1|stand,data=data)` The ANOVA for both models returns significant _p_ -values for moisture, but when I run the latter model, I get the following error message: `In pf(Fval[i], nDF[i], dDF[i]) : NaNs produced` ... the `NaN` is the _p_ -value for `stand`, which has zero denominator degrees of freedom.


title: Estimating demographics from internet activity

text: I'm looking for research papers on estimating demographic information (age, wealth, location, gender etc.) out of internet activities (preferably social network activity or browser history). Where would be a good place to look for?


title: Extracting lambda value with highest log-likelihood from boxcox output

text: library(MASS)     bc <- boxcox(Volume ~ log(Height) + log(Girth), data = trees)      To find the $\lambda$ value with the highest log-likelihood, this command could be used:               bc$x[which.max(bc$y)]      Is `x=log(Height)` in this example?


title: How to do post-hoc contrasts on within factor depending on certain between groups in SPSS?

text: Here are my experimental factors:   * within-factor (two levels): time1, time2.   * between factor (four levels): 4 groups. I want to test _specific post-hoc contrasts_ comparing:   1. group1: time1 vs. time2   2. group2: time1 vs. time2   3. group3: time1 vs. time2   4. group4: time1 vs. time2 to determine if and which group the time factor affects. SPSS _only_ provides an overall within effect for time averaging 4 groups. Moreover, SPSS _only_ provides post-hoc contrasts for the between groups averaging time How can I do my specific contrasts?


title: Predicting whether a potential sale will be won or lost

text: I am currently working on a project using a sales system and trying to come up with a way to use the current pipeline of potential sales to predict the amount of product that will be sold in the future. I’m looking for advice on how to approach this problem and hopefully some resources to teach me what approach to use and why. The sales system I’m using has historical data for opportunities (potential sales). Around 50,000 of the opportunities are “closed” meaning that they are either won or lost. I have around 1,000 “open” opportunities that have not yet been won or lost. Some variables that I have on each sale include the product (which is generally homogenous except for the amount), the amount, the salesman, the date, the time it was input into the system, the customer, and other data about the customer. I understand that if I want to predict a dichotomous variable like win / lose then I should look at a logistic regression. However, I’m looking for general advice on how to   1. Predict the probability of each individual opportunity closing as won using the data I have (and how to tell if I've done it correctly).   2. Estimate the total amount of won opportunities for a period. I found a similar question here Using a logistic model on the estimates of several other classification models but I’m hoping for a response that gives me a better idea of where to start. I’m comfortable using R or any other statistical software, but ideally I'd like some kind of book or other reference material that is as low-level as possible.


title: Quantifying differences in a discrete distribution across several populations

text: Suppose you have a discrete random variable, $Y$, with a large number (say, $300$) of discrete (which happen to be nominal) possible outcomes. The mass function, $p(y)=P(Y=y)$, is unknown but a sample is taken from $Y$ from 20 populations $i=1,2,...,20$. Let $Y_{ik}$ be the $k$'th sample from the $i$'th population and $n_{i}$ be the sample size taken from population $i$. Then we have the empirical mass functions: $$\hat{f}_{i}(y) = \frac{1}{n_i} \sum_{k=1}^{n_{i}} \mathcal{I}\{Y_{ik} = y\}$$ I'm trying to figure out a way of quantifying (and possibly testing) the difference in the distribution between populations. I've thought of one possible way to do this but I'd like to hear any other suggestions, even if they are completely different. What I've done so far is to create a measure for discrepancy between $\hat{f}_{i}$ and $\hat{f}_{j}$. One could treat this has a vector in 300-dimensional space and quantify differences by some metric (e.g. Euclidean distance). I wasn't too comfortable with this idea so, I create a matrix $$ D_{ij} = {\rm KS}(\hat{f}_{i}, \hat{f}_{j}) $$ where ${\rm KS}$ denote the Kolmogorov-Smirnov statistic. I then use $D$ as the distance matrix in multidimensional scaling, project the data down to 2 or 3 dimensions, and test heterogeneity between the points in this reduced dimensional space. This was a very ad-hoc approach and I'm guessing this type of problem has been worked on before so any advice or references would be much appreciated. **Note:** The sample sizes are small and vary between populations. The smallest populations may only have 5-10 samples (so there are MANY zero cells), while the larger ones may have more like 30.


title: How do I collect data for simple self-survey

text: I am on a self-experiment quest to solve some health related problems I have, and I want to create a self-survey that I take every day. What's the best way to codify the survey so I can later analyze the data? I am not a statistics guru, as you can tell, so I'm probably asking a pretty basic question here. For the survey itself, I plan to use a Google Docs Form. To give you some examples of the data I'd like to track, here are some example questions:   * Did I consume dairy?   * What was my stress level?   * How was my mood?   * How many hours of sunlight did I receive?   * Did I play with my children today? And so on.. My question is, how can I limit my own bias and quantify these things in a way that will be most useful?


title: What is the difference between EM and Gradient Ascent

text: I can't understand what is the difference between the algorithms EM (Expectation Maximization) and Gradient Ascent (or descent). Is there any condition under which they are equivalent? Thank you.


title: How to get dispersion parameter from a binomial mixed model?

text: I am modelling data with a generalized mixed model with binomial error distribution and I am concerned about overdispersion. I know that dispersion parameter can be measured as deviance/df, but for mixed models (= with random effect), the number of degrees of freedom cannot be extracted (I am using the function lmer from R). Is there any way to find out if my data are overdispersed?


title: classifiers providing probability of being correct

text: The way I understand classifiers like kNN, ANN, SVM, Decision Trees... is that after being trained they will associate a new test object to one class out of a set of predefined classes.      My question is whether there are classifiers, or extensions to the once mentioned above where objects are not only classified, but also a certain confidence level (probability of being correct) is provided. So something like $x $ belong to class $C_i$ with probability $0.8$.   I was especially hoping for something in connection with SVM, maybe based on the distance to the margin?


title: Drawing inferences from word counts

text: I'm trying to figure out if there is a good model that can connect word counts to classifying whether or not a person might go to an event. Here's the description: I have users, and events. 100 features in events are word counts; I don't know what each word is, but I know that feature 1 is the most used word stem of all event descriptions, feature 2 is the second most used word stem of all event descriptions... so on and so forth. I know which users go to which events. I just can't seem to find a way to incorporate these word count features into my current model. Does anyone know of a good direction to send me off in?


title: Correlated bivariate normal distribution: finding percentage of of data which is 2 standard deviations above the mean?

text: # Background Background: in an essay of mine, I point out that if a selection process (like higher education) requires successful applicants to be above the mean on 2 different variables, there will necessarily be fewer successful applicants than when successful applicants just had to be above the mean on just 1 variable. (The 2 different variables in this case are the Big Five personality factor of Conscientiousness, and IQ.) This reduction is most dramatic when the 2 variables have 0 correlation, but even a large correlation will still result in many applicants being filtered out. How many are filtered out, exactly? # Simple questions Well, if the filter is for 2 standard deviations above the mean and the variables are correlated with 1, then 2.3% of the population will pass; if the variables are uncorrelated with 0, then 2.3% * 2.3% or 5.29e-4% of the population will pass. # Correlated But what about intermediate values? For example, the psychology literature has reported a correlation of -0.21 between Conscientiousness & IQ. I consulted Wikipedia on bivariate normal distributions, but I didn't understand much of it. The closest I found was sum of correlated normal random variables, but in this case what I want is closer to a `min` function. ## Simulation I wa able to work up a R simulation to see how that worked, and it seemed in line with my intuitions:               install.packages("fMultivar")     library ("fMultivar")          x <- rnorm2d(10000000, rho=0.5)          xgreater <- length(subset(x, x[,1] > mean(x[,1])+2*sd(x[,1])))     xandygreater <- length(subset(x, x[,1] > mean(x[,1])+2*sd(x[,1]) & x[,2] > mean(x[,2])+2*sd(x[,2])))          c(xgreater, xandygreater); c(xgreater / length(x), xandygreater / length(x), xgreater / xandygreater) * 100          # example results for different values of 'rho='     0.1     [1] 454,664  17,570     [1] 2.273e+00 8.785e-02 2.588e+03          0.2     [1] 458,284  82,552     [1]   2.2914   0.4128 555.1458          0.5     [1] 454,484  80,872     [1]   2.2724   0.4044 561.9794          0.9     [1] 455,242 267,912     [1]   2.276   1.340 169.922          0.95     [1] 455,162 321,024     [1]   2.276   1.605 141.784          0.99     [1] 455,260 394,448     [1]   2.276   1.972 115.417      ## Exact pdf calculation? I really was hoping for more of a precise analytic solution, so some more searching eventually turned up a paper, "Exact Distribution of the Max/Min of Two Gaussian Random Variables", which gives a definition for the min of 2 correlated normal variables. This seems to be what I want; top of pg1, second column: > ...where $\phi(.)$ and $\Phi(.)$ are, respectively, the pdf and the > cumulative distribution function (cdf) of the standard normal distribution. > It is known that the pdf of $Y = \min(X_1, X_2)$ is $f(y) = f_1(y) + > f_2(y)$, where > > (3) $f_1(y) = \frac{1}{\sigma_1} \phi (\frac{y-\mu_1}{\sigma_1}) \times \Phi > (\frac{p(y - \mu_1)}{\sigma_1 \sqrt{1 - p^2}} - \frac{y - \mu_2}{\sigma_2 > \sqrt{1 - p^2}})$ (4) $f_2(y) = \frac{1}{\sigma_2} \phi > (\frac{y-\mu_2}{\sigma_2}) \times \Phi (\frac{p(y - \mu_2)}{\sigma_2 \sqrt{1 > - p^2}} - \frac{y - \mu_1}{\sigma_1 \sqrt{1 - p^2}})$ They give an R implementation on pg6 (first column); it seems to have a `pnorm` typo, but I fixed that. Once it was working, I tried generating a slightly (0.1) correlated bivariate distribution, which look OK:               fmin<-function (y,mu1,mu2,sigma1,sigma2,rho)          {t1<-dnorm(y,mean=mu1,sd=sigma1)          tt<-rho*(y-mu1)/(sigma1*sqrt(1-rho*rho))          tt<-tt-(y-mu2)/(sigma2*sqrt(1-rho*rho))          t1<-t1*pnorm(tt)          t2<-dnorm(y,mean=mu2,sd=sigma2)          tt<-rho*(y-mu2)/(sigma2*sqrt(1-rho*rho))          tt<-tt-(y-mu1)/(sigma1*sqrt(1-rho*rho))          t2<-t2*pnorm(tt)          return(t1+t2)}     fmin(c(1:200),100,100,15,15,0.1)       [1] 1.849e-11 2.864e-11 4.418e-11 6.784e-11 1.037e-10 1.578e-10 2.392e-10 3.608e-10 5.418e-10      [10] 8.101e-10 1.206e-09 1.787e-09 2.636e-09 3.872e-09 5.663e-09 8.243e-09 1.195e-08 1.724e-08      [19] 2.476e-08 3.542e-08 5.043e-08 7.148e-08 1.009e-07 1.417e-07 1.982e-07 2.760e-07 3.827e-07      [28] 5.282e-07 7.257e-07 9.928e-07 1.352e-06 1.833e-06 2.475e-06 3.326e-06 4.449e-06 5.926e-06      [37] 7.859e-06 1.037e-05 1.364e-05 1.784e-05 2.324e-05 3.014e-05 3.891e-05 5.002e-05 6.401e-05      [46] 8.154e-05 1.034e-04 1.306e-04 1.641e-04 2.054e-04 2.558e-04 3.173e-04 3.917e-04 4.814e-04      [55] 5.889e-04 7.173e-04 8.696e-04 1.049e-03 1.261e-03 1.507e-03 1.794e-03 2.125e-03 2.506e-03      [64] 2.941e-03 3.435e-03 3.993e-03 4.620e-03 5.318e-03 6.093e-03 6.945e-03 7.878e-03 8.890e-03      [73] 9.982e-03 1.115e-02 1.239e-02 1.370e-02 1.506e-02 1.647e-02 1.791e-02 1.938e-02 2.084e-02      [82] 2.230e-02 2.371e-02 2.508e-02 2.636e-02 2.755e-02 2.863e-02 2.956e-02 3.034e-02 3.095e-02      [91] 3.138e-02 3.162e-02 3.165e-02 3.149e-02 3.112e-02 3.056e-02 2.981e-02 2.889e-02 2.781e-02     [100] 2.660e-02 2.526e-02 2.383e-02 2.233e-02 2.077e-02 1.920e-02 1.762e-02 1.605e-02 1.452e-02     [109] 1.305e-02 1.164e-02 1.031e-02 9.063e-03 7.912e-03 6.857e-03 5.899e-03 5.039e-03 4.272e-03     [118] 3.595e-03 3.004e-03 2.491e-03 2.050e-03 1.675e-03 1.358e-03 1.093e-03 8.732e-04 6.923e-04     [127] 5.447e-04 4.254e-04 3.297e-04 2.535e-04 1.935e-04 1.466e-04 1.102e-04 8.220e-05 6.085e-05     [136] 4.470e-05 3.258e-05 2.357e-05 1.692e-05 1.205e-05 8.517e-06 5.973e-06 4.157e-06 2.870e-06     [145] 1.966e-06 1.337e-06 9.018e-07 6.036e-07 4.008e-07 2.641e-07 1.727e-07 1.120e-07 7.210e-08     [154] 4.604e-08 2.917e-08 1.834e-08 1.144e-08 7.078e-09 4.346e-09 2.647e-09 1.600e-09 9.594e-10     [163] 5.708e-10 3.369e-10 1.973e-10 1.146e-10 6.607e-11 3.778e-11 2.143e-11 1.207e-11 6.737e-12     [172] 3.733e-12 2.052e-12 1.119e-12 6.052e-13 3.248e-13 1.730e-13 9.137e-14 4.788e-14 2.489e-14     [181] 1.284e-14 6.571e-15 3.336e-15 1.680e-15 8.393e-16 4.160e-16 2.046e-16 9.979e-17 4.830e-17     [190] 2.319e-17 1.104e-17 5.218e-18 2.446e-18 1.137e-18 5.248e-19 2.402e-19 1.090e-19 4.911e-20     [199] 2.194e-20 9.727e-21      ![Plotting the pdf](http://i.imgur.com/10n1FuI.png) Now, I understand the PDF to be "a function that describes the relative likelihood for this random variable to take on a given value. The probability for the random variable to fall within a particular region is given by the integral of this variable’s density over the region". So I suppose I should sum up every point in the pdf >130 (since 130 is 2 standard deviations up, by construction when I specified SD=15) and that's my probability that a random deviate will be `min(130,130)`. What's the total probability someone will be over 130 on both variables? I think that would be:               sum(fmin(c(1:200),100,100,15,15,0.1)[130:200])     [1] 0.001004      If I increase the _r_ to 0.9, the result is 0.01455 which is satisfyingly larger. A sanity check - as the correlation goes to 1.0, there should be no decrease. So we do the same question for a single normal distribution defined the same way:               sum(dnorm(c(1:200), 100, 15)[130:200])     [1] 0.02459          # the function blows NaN chunks on 1.0, so we'll try a lot of 9s:     sum(fmin(c(1:200),100,100,15,15,0.9999999999)[130:200])     [1] 0.02459      I guess that works too. # Problems So my questions are:   1. Is my R simulation right?   2. Is my version and use of `fmin` right?   3. Is there some more direct, possibly even pen-and-paper, avenue of calculating the answer to my original question?


title: Likelihood ratio test for sex ratio

text: I am trying to solve this problem: Assume that each child is male with probability $p$ independently of all other children. We observed 19711 male births out of a total of 38562 births in American families with two children each. Use the likelihood ratio test statistic to test the hypothesis $H_0: p=1/2$ against a suitable alternative, which you should specify. We also found 17703 males out of 35042 similar births in Finland. Use the generalised likelihood ratio test to test the hypothesis that $p$ has the same value in each country versus a suitable alternative. For the first part we have $X_1,..X_{38562}$ I.I.D. r.v. with distribution $\mathrm{Ber}(p)$. The alternative hypothesis should be $H_1 : p \in [0,1]$. Now as we know that the sample mean is the MLE for Bernoulli distributions, we have likelihood ratio, with $n= 38562$: $$ \displaystyle \lambda(X)= \frac{(1/2)^n}{\overline x^{\sum x_i} \cdot (1-\overline x)^{(n-\sum x_i)} } $$ So the likelihood ratio statistic is: $$ \Lambda(X) = -2 \log(\lambda(X)) $$ Hence our test is that we reject $H_0$ iff $\Lambda(x) \geq k$ for a constant $k$ such that $$ P(\chi^2_1 \geq k) = \alpha. $$ Where $\alpha$ is the required size. We can do that as $n$ is large and we know that the likelihood ratio statistic converges in distribution to the chi squared with 1 degree of freedom ( **QUESTION** : the degrees of freedom should be $\mathrm{Dim}([0,1])- \mathrm{Dim}({1/2})= \mathrm{Dim}([0,1])$ why do we assign dimension 1 to this? I think we do, because otherwise 0 degrees of freedom wouldn't make much sense...) I thought a sensible choice of $\alpha$ would be 0.05 so I worked out $k= 3.841$. How do I proceed now? Putting the values in the likelihood ratio statistic leads to impossible calculations... Any advice on this last bit or on the following part would be really appreciated!!


title: Repeated measures ANOVA, with same subjects under 2 treatments

text: I have a group of 6 animals. I took the baseline reading (Time 0), put them under Treatment A, then take readings again at Times 2, 4, 6, and 12 hours. I rest them for 2 weeks. Then I took a baseline (time 0), put them under Treatment B, then take readings again at Times 2, 4, 6, 12 hours. I believe this is a repeated measures ANOVA, with within-subject (as Times- with 5 levels: 0,2,4,6,12 hrs) and between-subject (Treatment A, Treatment B). My problem is, the subjects in Treatment A and Treatment B are the same subjects. How do I account for that (as in they are not randomized)? Right now, I ran the repeated measures MANOVA and ANOVA for time*treatment, time, treatment. Then I ran paired-t test with bonferroni correction (0.05/4) to test Times 2,4,6,12 with Time 0 (of the same treatment). And I ran paired-t test to test Treatments (A vs B) across time points (Time 0,2,4,6,12).


title: R: statistical test to identify samples with too high variability

text: I would like to develop a test to identify which variables in my data set have a variation higher than the "average variability". I'm struggling with that since days, and I also tried in vain to look for help in other forums. I have data from biological experiments, that look like this:               v1 2 1.8 1.5 1.9 2.1 1.78 1.95 2.0 2.1       v2 2 100 -5.2       v3 1 -1.3 -2 2.3       v4 1 1.5 1.6 1.9 2.1 2.0 2.4 -1.1 2.3 1.5 1.6 1.9 1.8 1.6      These represent gene expressions. Now, I would expect that all values of each variable(genes) are more or less similar, since the values are ripetute measurements of the same gene. Having a variable with such a huge difference, as v2 , doesn't have sense, because the repeated measurements should give consistent values. Therefore, it has to come from a methodological error and the variable(gene) has to be discarded. I was looking for a method (possible a statistical test) in R which could identify the "average variability" among my samples and report me which variables(genes) have a variability significantly greater. This means that for these genes my data are not good enough to estimate the expression, and I have to discard them. I would really appreciate any suggestion/links/advice/methods on test I could use for my purpose.


title: Covariance between f(x+y) and f(x)

text: Is there some theorem that allows us to work out: $Cov\big( f(x+y), f(x)\big)$ We now the $Var\big(f(x)\big)$ and $Var\big(f(x+y)\big)$ and also the $Var(x)$, $Var(y)$ and the $Covar(x,y)$ Our specific case is the following: What is the $Covar\big(\Phi(\beta_0+\beta_1),\Phi(\beta_0 )\big)$, with $\beta_0$ and $\beta_1$ probit model coefficients with variance and covariance known.


title: Derivation of the Simpson index

text: I'm interested in the derivation of the Simpson diversity index. I found the normal interpretation of this index, but I'm not able to see what's the "mathematical background" of the formula. I read the original paper of Simpson, but I'm more interested in an easier explanation (my skills in probability theory are more on an undergraduate level). I would appreciate a little help here. Thanks in advance Meiner


title: Obtain the coordinates of one point (X1,Y1) where it is maximum the change of the curve

text: ![a busy cat](http://www.biomedcentral.com/content/figures/1471-2199-9-66-1-l.jpg) Dear community, I used a figure obtained from internet to illustrate what I would like to solve. Imagine this exponential regression curve. How can I get the exact point of the curve where it is maximum the change of the it? Observing the figure, it should be placed between 1 and 2 in the X axis and between 0.04 and 0.08 in the Y axis. But where it is exactly? Thanks in advance


title: Should I use a beta-binomial or binomial glmm?

text: I have several data sets on wildlife disease incidence. One of the issues with my dependent variable is that it represents only current infection status, therefore 0 (no disease) can represent either an animal who has never had the disease or has had it and recovered (it is an acute viral infection). There is evidence of overdispersion in my datasets. Therefore I after having already performed model selection on a binomial glmm (in R, `lme4` package), I decided that perhaps beta-binomial model would be preferable. I fitted all subsets of models in R's `hglm` package and performed model selection based on AICc (I believe this is calculated from the h-likelihood). The results of this model selection process are quite different to those from the binomial glmm, as one would expect, however I'm inclined to be a bit suspicious of the results from the beta-bin model, as in all 4 cases the global model was chosen as best - so perhaps this measure is just not parsimonious enough? Am I using the right information criteria?


title: What laws direct marketing obeys?

text: I am starting a new product. I need 100 clients to be viable, out of a market of 12000 target companies. I'm doing direct marketing in 2 steps: Call / Emailing. After 50 calls, I collected 30 emails. Out of 30 emails, I got 37.5% opens, 13.6% clicks, and 0% new customers. But where should I stop? My guess is the following: I need 1% of the market. I should expect 1% returns. This calls for 100 attempts leading to 1 return. But I'm lost right here: What is the accuracy of 100 attempts? Does it say: "You have 0 returns on 100 attempts, you underdo your target with 90% certainty", or "95% certainty", or "50% certainty"? I feel I need at least 90% certainty. So, where should I stop? Thank you


title: P-value of F-test to compare two variances (var.test in R)

text: I am trying to understand where the p-value of a F-test comparing two variances comes from. More specifically, the p-value given by R's `var.test` function does not match p-value assigned to a F-test by the `pf` function with the same F value and degrees of freedom. For example, p-value given here:               > d1 <- rnorm(300, sd=1)     > d2 <- rnorm(300, sd=1.2)     > var.test(d1, d2)              F test to compare two variances          data:  d1 and d2     F = 0.78, num df = 299, denom df = 299, p-value = 0.03212     alternative hypothesis: true ratio of variances is not equal to 1     95 percent confidence interval:      0.62 0.98     sample estimates:     ratio of variances                    0.78       Does not match this one:               > pf(0.78, 299, 299, lower.tail=F)     [1] 0.98      Could someone explain where the difference comes from?


title: Outliers in a bimodal distribution

text: _Somewhat related tothis post._ However my question is a little specific. I have some data from gene expression quantification from single cells using RNA sequencing. I haven't had a look at it but I suspect that there may be errors in the procedure that may lead to high variability. Also, it is quite possible that certain genes may exhibit multimodal (at least bimodal) expression owing to stochastic fluctuations or bistability. In such cases, would this be a right way of removing outlier(s):   * Generate distribution    * Remove bins (with data) which have a probability less than a threshold $\alpha$ I am also unsure of whether there is any logical (formulaic) way of deciding bin sizes.


title: Recursive logistic regression merge

text: I need to make regression on big amount of data, each row have around 1000 features. Did will outcome will be same or better when i make 4 separate regressions of 250 features and after that i will make one regression that will have 4 features equal to underlying regression outputs? I can't make one regression on all features because coefficient learning algorithm using to much memory. General known solution: $ Y' = sigmoid(X*\beta) $ My take: $ X' = $\begin{bmatrix} Y_1' = sigmoid(X_{1-250}*\beta_1) \\\ Y_2' = sigmoid(X_{251-500}*\beta_2) \\ Y_3' = sigmoid(X_{501-750}*\beta_3) \\ Y_4' = sigmoid(X_{751-1000}*\beta_4)\end{bmatrix} $ Y'' = sigmoid(X'*\beta) $ I'm asking about differences/relationships between $Y'$ and $Y''$. Rows (observations) around 100 000 000, features 1000 ($Length[X]$) Sorry for bad formatting in equations, I'm not a MathJax master.


title: Does the importance of the normality assumption in MANOVA decrease as sample size increases?

text: Does the importance of the normality assumption in MANOVA decrease as group sample size increases?


title: Is the term "Shannon filter" appropriate for a naive filter based on the Shannon information entropy of a data set?

text: This relatively simple question collapses neatly into the title. I have large swaths of numerical and categorical data that contain unknown amounts of signal relative to weak data modeling principles. From these data sets, I'm attempting to extract usable signal -- that is, data that is well-ordered, correctly transformed and modelled, and useable for stochastic ensemble construction and analysis that occurs downstream. In order to do this, I have decided to construct what I am calling a _Shannon Filter:_ a simple, naive system that computes the Shannon information entropy of each piece of data with respect to its current model, ranks the data in descending order (relative to a logitnormal distribution), and performs filtering at some confidence level yet to be determined. While this is not yet a methods question (I may turn to http://dsp.stackexchange.com or http://crypto.stackexchange.com for that answer, since both spaces contain direct analogs to my problem), it is currently one of correct terminology. The phrase "Shannon Filter" appears to be in loose current usage, yet I can find no sourcing to back this up. Would use of this term be considered appropriate in the statistical domain? * * * _I ask here, because of all the StackExchange sites that delve into this class of problem (inclusive of math, dsp, crypto, cs-theory, SO, and english; yikes!), abstract statistics seems to be the most relevant domain that I wish to apply this term within. I expect my system to be defined in the context of stochastic sampling and pairwise search algorithms digestible by statisticians, so I require a crisp, clear definition for this process that is technically correct._ _Thank you in advance._


title: I do not understand the terminology of $Q_n$

text: I am trying to understand what looks to be some simple terminology issue, involving $Q_n$ from my book here: ![enter image description here](http://i.stack.imgur.com/MiVuH.jpg) Specifically, I do not understand what they mean when they say $m = {n \choose 2}$, and what they mean when they say $k = {[n/2] + 1 \choose 2}$. As far as I know, this terminology means "something choose something", (in relation to combinations, etc), but I do not see how it possibly applies in this case. Am I misinterpreting this terminology, or does it mean something all together, and if so, what does it mean? Thanks in advance.


How can I calculate proportion of variance of an aggregate time-series accounted for by the variance of a sub-series? As I'm not sure what the correct terminology is for explaining my question, I'll illustrate it with an example: I have a time-series (say, for example, annual alcohol consumption in the USA) that is exactly made up of the sums of several sub-series (continuing the example, annual consumption of spirits, beer and wine). The variation in one of these component time-series is much greater than in the others, meaning that it dominates the variation in the aggregate series. For example:               n.obs <- 20     set.seed(1990)          x1 <- cumsum(rnorm(n.obs, mean = 0, sd = 1)) + 20     x2 <- cumsum(rnorm(n.obs, mean = 0, sd = 1)) + 20     x3 <- cumsum(rnorm(n.obs, mean = 0, sd = 1)) + 20     x4 <- cumsum(rnorm(n.obs, mean = 0, sd = 8)) + 100          y <- x1 + x2 + x3 + x4      As you can see, the variation in `y` is dominated by the variation in `x4`: ![Plot of y, x1, x2, x3 and x4 series overlayed. y visually follows the path of x4, but with a higher intercept](http://i.stack.imgur.com/s4bKW.png) **How can I quantify the extent to which the variation of`y` is accounted for by the variation in `x4`?** My initial thought was to calculate the variance in `y-x4` compared to the variance in `y` (i.e. `1 - (var(y - x4) / var(y))`) but if I do this for all sub-series the values do not sum to one as I would expect.


title: How do I transform a data generation specification in R into a BUGS/JAGS specification

text: I am at a loss about what the BUGS/JAGS specification of the following should look like. The background is that one person has four measurements taken by four different instruments. Each instrument is expected to have a different associated error. The R code for generating simulated data (not very elegant) is:               n <- 63                    # Number of patients          ## (Hyper)parameters     var_a <- 2                  # Variance around population mean     mu_a <- 80                  # population mean     mu_m <- c(0, -5, 1, -6)     # Mean diff for each measurement instrument from the population mean     var_m <- c(1, 1.5, 2, 2.5)  # Variance of each measurement instrument     m_i <- 0                    # Individual measurement error (mean)     m_e <- 1                    # Individual measurement error (variance)          # Counter to keep track of total number of measurements taken     measures.taken <- 0               ## Generate the data     for (i in 1:n){                                                       # loop around the number of patients       a <- mu_a + rnorm(1, 0, sqrt(var_a))                                # Random intercept       for (j in 1:length(mu_m)){                                          # Loop to generate 4 measures for each patient         measures.taken <- measures.taken + 1                              # Just keeping track of the number of measurements in total         individual.error <- rnorm(1, m_i, sqrt(m_e))            bp <- rnorm(1, mu_m[j], sqrt(var_m[j])) + a + individual.error   # This sets up the random variability for the measures         print(c(bp, i, j))                                               # Print data as generated to the standard output       }      }      Any insights gratefully accepted.


title: Why is it bad to teach students that p-values are the probability that findings are due to chance?

text: Can someone please offer a nice succinct explanation why it is not a good idea to teach students that a p-value is the prob(their findings are due to [random] chance). My understanding is that a p-value is the prob(getting more extreme data | null hypothesis is true). My real interest is **what is the harm** of telling them it is the former (aside from the fact it is simply not so).


title: Methods to check if my data fits a distribution function?

text: My shortened data is:               y <- c (2,2,1,5,6,7,1,2,1,6,6,7,3,2,4,4,4,4,3,3,9,1,1,9)      I firstly normalize my data:               y_scale <- scale(y)      Then, I generate a model dataset with normal distribution based on y_scale's mean and stdev:               y_norm <- rnorm(n=24, m=mean(y_scale), sd=sd(y_scale))      To check if my data fits the normal distribution, I do               ks.test(y_scale,y_norm)      I found the result is as follows:               Two-sample Kolmogorov-Smirnov test          data:  y_scale and y_norm      D = 0.2083, p-value = 0.6749     alternative hypothesis: two-sided           Warning message:     In ks.test(y_scale, y_norm) : cannot compute correct p-values with ties      Here, my question is: (1) My real data set has ~ 700,000 numbers, I found I cannot use shapiro.test.               shapiro.test(y_scale)     Error in shapiro.test(y_scale) : sample size must be between 3 and 5000      (2) Is the p-value calculated above by ks.test wrong? How to solve this problem of p-values?               Warning message:     In ks.test(y_scale, y_norm) : cannot compute correct p-values with ties      (3) The reasons why I tried to use ks.test instead of other methods, is because I want to compare with other model datasets that have other distribution functions. It seems to me that I can simply replace y_norm with other model dataset, and compare their p-values or D-values (the smaller the better), to choose which distribution function fits my data most. (4) Is it a must to normalize my data first?


title: Regular data gaps in a time series

text: I've been looking at attendance data for a college library taken during finals week when the library was open for 24 hours. However the staff only took head counts during the hours of 12am to 6am (7 data points per night, for about 11-12 nights, 3 semesters worth of similar counts). I was wondering if, because the data isn't "continuous" I can really make anything analysis on it, since it seems to create a cyclical pattern different from what I would see over a true 24 hour period. Do you think the data is still usable? Also if anyone knows of any special tool applicable to time series analysis using attendance data, where we want to know how many are using the library and how to get more people to attend, I would be extremely grateful, as I'm new to statistical analysis in general. Thanks!


title: How to tune smoothing in mgcv GAM model

text: I am trying to figure out how to control the smoothing parameters in an mgcv:gam model. I have a binomial variable I am trying to model as primarily a function of x and y coordinates on a fixed grid, plus some other variables with more minor influences. In the past I have constructed a reasonably good local regression model using package locfit and just the (x,y) values. However, I want to try incorporating the other variables into the model, and it looked like generalized additive models (GAM) were a good possibility. After looking at packages gam and mgcv, both of which have a GAM function, I opted for the latter since a number of comments in mailing list threads seem to recommend it. One downside is that it doesn't seem to support a local regression smoother like loess or locfit. To start, I just wanted to try to replicate approximately the locfit model, using just (x,y) coordinates. I tried with both regular and tensor product smooths:               my.gam.te <- gam(z ~ te(x, y), family=binomial(logit), data=my.data, scale = -1)          my.gam.s  <- gam(z ~  s(x, y), family=binomial(logit), data=my.data, scale = -1)      However, plotting the predictions from the model, they are much much more smoothed compared to the locfit model. So I've been trying to tune the model to not oversmooth as much. I've tried adjusting the parameters sp and k, but it's not clear to me how they affect the smoothing. In locfit, the nn parameter controls the span of the neighborhood used, with smaller values allowing for less smoothing and more "wiggling", which helps to capture some areas on the grid where the probability of the binomial outcomes changes rapidly. How would I go about setting up the gam model to enable it to behave similarly?


title: Z-scores, standardized tests and population means

text: I have been using z-scores by subtracting the mean (and then dividing by the SD) of the _sample_ , whereas I recently read they actually need to be the _population_ mean and SD. Assuming you are transforming raw scores to z-scores in order to compare within-sample variables (various test scores of the same sample of subjects), does it really invalidate the comparison if you have been using sample rather than population statistics? Also, is it _ever_ realistic that one would know the population statistics? I always assumed this is more of a theoretical concept than an actual statistic. For instance, in order for the population mean of the IQ test to be set to 100, surely this didn't imply testing every single human on the planet and adjusting the raw-to-standardized tables to reflect the population average! Many thanks.


title: Linear regression: Evaluate probability of $Y>y| X=x$

text: Given a linear regression model with all the assumptions checked and validated, I would like to obtain the probability that $Y>y|X=x$. For example for the iris dataset, I would do the following to obtain the probability of $Y>5|X=1,2,3...7$:               plot(Sepal.Length~Petal.Length, data=iris)     lm1<-lm(Sepal.Length~Petal.Length, data=iris)     summary(lm1)     abline(lm1)     predict(lm1, newdata=data.frame(Petal.Length=1:7))     (summary(lm1))$sigma         pnorm(5, mean = predict(lm1, newdata=data.frame(Petal.Length=1:7)),             sd = (summary(lm1))$sigma, lower.tail = F)      Is such an approach correct assuming constant variance?


title: neural net reinforcement learning board game, how to achieve higher competence?

text: I'm trying to write an algo for playing some board games using Q/SARSA combined with neural net as valuation approximation. It works to some extent, but I can't make it to expert level. So for example, my AI can play tic-tac-toe perfectly, or the game 2048 perfectly on 3x3 board. But then if I ask it to learn connect-4 or the full 2048, it can only learn to a very basic level. Eg play 2048 to 128/256, but not higher. I guess my implementation is correct as it learns easier games without problem, but as the game becomes complicated, the learning stopped at basic level. Can someone point out references to go further? Or any suggestions? Thanks. Implementation detail: I followed Sutton's TD-gammon example. Eg connect 4, I have a 1-hidden layer feed forward net, with 128 input neurons (64 squares, each square takes 2 neurons.) I tried 40/80/120 neurons for hidden. Hidden/Output layers with tanh activation. SARSA learning rule. External reward is 0 unless game finished. I make it play against itself for >1million games.


title: Can I normalize a data proportion?

text: Is there any way to normalize this variable (attached)? In fact, this is a proportional data (a percentage). But I need transform it in order to do some contrasts. I try with arcsine transformation: 2 * ARSIN(SQRT(x)). But I can't get an interesting result. Moreover, If there are any way to get a normal distribution from this data, how should I report my results? For mean comparisions, can I report observed differences in a non transformed variable and accompany with the p-values derived from the transformed variable?![enter image description here](http://i.stack.imgur.com/wr0Yj.png)


title: Probability of 3 brothers born on the same date of 3 different months?

text: Whats the probability of 3 brothers being born on the same day of the month in 3 different months? I ask because I'm one of them. Its April, May, and October 12th. What's the difference including whether or not we'd be male/female?


title: Metric to show the overall similarity between two set of points obtained from handwritten word?

text: In my case, I would like to compare different handwriting to a template handwritten word, and return a single metric to show how closely the input matches the template's shape. The lines are the same thickness; their main difference should be the curve and angle. I have already collected and paired up two sets of points from two line images. Then my next step is to sum up a single metric to describe the similarity/variance in these 2D images. I learnt about Cosine similarity, Euclidean Distance, Minkowski distance, variance of distance, etc. However, I am not sure of their limitations in my case (may lower accuracy). Please suggest which one should I go for. Thanks.


title: Recovering true data from many noisy samples with varying unknown amounts of noise

text: **Input:** $k$ vectors $x^1,\ldots,x^k \in \mathbb{R}^n$, where $x^i \sim \mathcal{N}(x,\mathbb{1} \cdot \sigma_i^2)$. **Goal:** approximate the vector $x$ as well as possible. The quality of approximation is measured in any reasonable metric such as $\mathbb{L}_1$ norm or Pearson correlation. Note that the values of $\sigma_i$ are not given as part of the input. In my input, $k \gg n$. What is a good solution to this problem? I'm pretty sure one can do better than simply returning the average of all the vectors $x^i$. One possibility I'm thinking of is performing a variant of canonical correlation analysis between $x^1$ and $x^2,\ldots,x^k$: the variant I'm thinking of would find the best linear combination only among those that have _only positive coefficients_. (I think this can be solved using semidefinite programming.)


title: Estimating Negative Binomial Regression Model

text: It is easy to estimate a Poisson regression model using the Newton–Raphson Iterative Technique as it only involves one parameter (mu). However, I am unable to understand how a negative binomial regression model can be estimated using the Newton–Raphson Iterative Technique, as it involve two parameters (mu and dispersion parameter alpha). Can anyone help me in this regard?


title: How can we evaluate the predicted values in regression model

text: I am using regression model to predict values I have. How can evaluate the accuracy of prediction model (I'd like to see how the accuracy of predicted values). I found different metrics, but it's only used with discrete data.   * Classification report   * Accuracy score   * Confusion matrix Now, I can generate the summary using .summary() function. As an example here , it's nearly the same in R/Python; I'd like to evaluate the performance, how good and accurate the predicted values that I got using GLM (with Binomial families ).


title: How do multi-attribute edge-weights influence community detection?

text: My graph consists of a computer network topology where each vertex is a physical node/device (depicted using its IP address). Two vertices will have an edge if the nodes have had communication with each other. We call each communication a _conversation_ , which will have multiple packets transferred in both directions. Multiple edges between two nodes are thus obvious. For any _conversation_ , I need to capture multiple features, such as:   1. Bytes transferred between the two nodes   2. No. of packets transferred between the nodes   3. The duration of data transfer between the nodes   4. The inter-arrival time of packets between the nodes While trying to find communities in my graph, I am coming across certain difficulties:   1. Should I go for multi-attribute edges (with each of the above feature taken as a weight), or should I be taking something like weighted average of all features and make them a single attribute?   2. Will I obtain better results if each of the above feature is taken as a separate edge-attribute? If so, how and why?. Does it influence 'modularity' calculations? Again, how? (Links to past works or papers will do)   3. What algorithms may suit my task well? I do not know about other algorithms apart from Louvain (igraph implementation).


title: Help required with using the correct model for sampling

text: I regularly encounter the following scenario: I am presented with a set, $\mathcal{S}$, of $i$ items. I am assured that $d$ of the $i$ items possess a particular property, $d'$ - and that the rest, $(i-d)$, do not. So the items are in 2 piles: $d'$ and non-$d'$. (For context, let's say that $i$ is known to be 1,000,000 and $d$ is said to be 30,000. ) I stipulate that the probability of any one member of $i$ also being a member of $d$ is independent of other members of $i$. Consider the case where I wish to test this assurance by taking a sample. (For practical and financial reasons, I cannot test all of $\mathcal{S}$. ) Clearly I don't know how common the property of $d'$ actually is within $\mathcal{S}$.   1. What size of sample from $\mathcal{S}$ do I need, for a given confidence and margin of error to test whether the provided value of $d$ is correct or not? E.g., if I was told $d$ was 30,000 and my sample indicated it was 30,100 I would be happy. If the indication was that the real value was 40,000 I would not be happy.   2. Suppose I only have access to the items marked $d'$ and cannot sample the rest. I believe then I can only test for non-$d'$ items in the pile marked $d'$, what I call false positives. Can I then infer anything about false negatives, e.g. $d'$ items in the non-$d'$ pile which should not be there?   3. Suppose further that in fact some proportion of the items in $\mathcal{S}$ do have a dependency, such that if a particular item $p$ has $d'$ then item $p+1$ has $d'$ also . Could I reasonably remove the dependency by considering i and i+1 as a single item? If someone could point me at the relevant area of statistics, I would be very grateful!


title: Intuitive explanation of Fisher Information and Cramer-Rao bound

text: I am not comfortable with Fisher information, what it measures and how and how is it helpful. Also it's relationship with the Cramer-Rao bound is not apparent to me. Can someone please give an intuitive explanation of these concepts?


title: Testing significant effect in 2 by 2 factor design on a binary outcome

text: I have 4 groups each consisting of 1000 recent college graduates, and I am studying employment immediately after graduation. 2 groups took a business class in college, and 2 groups took a math course in college. Group 1 is a control that has taken neither a business class nor a math course in college. I want to conduct a statistical analysis to find if either class was significant in increasing the graduate in obtaining a job immediately after school. Here's the data: Group 1: No Math, No Business 1000 subjects, 50 employed Group 2: No Math, Yes Business 1000 subjects, 60 employed Group 3: Yes Math, No Business 1000 subjects, 75 employed Group 4: Yes Math, Yes Business 1000 subjects, 100 employed


title: How do I propagate error values through a matrix diagonalization?

text: In the general case, if I have a symmetric, invertible (real-valued) matrix where each entry in the matrix has a separate error associated with it (or, in the more general case, where each value has a separate positive and negative error associated with it), how can I _rigorously_ propagate these error values through the eigendecomposition to get error values for each of the eigenvalues and each entry in the eigenvector matrix? For my current application, it is as critical that I avoid overestimating the error as it is that I avoid underestimating it.


title: Power calculation for likelihood ratio test

text: I have two independent poisson random variables, $X_1$ and $X_2$, with $X_1 \sim \text{Pois}(\lambda_1)$ and $X_2 \sim \text{Pois}(\lambda_2)$. I want to test $H_0:\, \lambda_1 = \lambda_2$ versus the alternative $H_1:\, \lambda_1 \neq \lambda_2$. I already derived maximum likelihood estimates under null and alternate hypothesis (model), and based on those I calculated likelihood ratio test (LRT) statistic (R codes given below). Now I am interested to calculate power of the test based on:   1. Fixed alpha (type 1 error) = 0.05.   2. Using different sample sizes (n), say n = 5, 10, 20, 50, 100.   3. Different combination of $\lambda_1$ and $\lambda_2$, which will change LRT statistics (computed as `LRTstat` below). Here is my R code:               X1 = rpois(λ1); X2 = rpois(λ2)     Xbar = (X1+X2)/2     LLRNum = dpois(X1, X1) * dpois(X2, X2)     LLRDenom = dpois(X1, Xbar) * dpois(X2, Xbar)     LRTstat = 2*log(LLRNum/LLRDenom)      From here, how could I proceed with power calculation (preferably in R)?


title: Why do we refer to our estimates in terms of precision?

text: Open any statistics textbook and it will urge the need to check the 'precision of our estimates'. Take the following random variable:               set.seed(55)     x <- rnorm(n=10, mean=0.3453, sd=1)          x      [1]  0.4654391 -1.4670769  0.4968830 -0.7739210  0.3472082  1.5338185 -0.1600439  0.2460656  0.6506532     [10]  0.5437097          mean(x)     [1] 0.1882736          t.test(x ,alternative="two.sided",mu = 0,conf.level=0.95)$conf.int     [1] -0.4025211  0.7790682      I would describe the mean of this variable as being very precise, in fact precise to 7 decimal places. I would not describe it as being very accurate, as its 95% confidence interval spans a large range of the data. However, in statistics-speak, the mean of the variable would be described as being imprecise. So why don't we refer to the 'accuracy of our estimates', instead of the 'precision of our estimates'?


title: Bayesian analysis with histogram prior. Why draw simulations from the posterior?

text: This is a beginner’s question on an exercise in Jim Albert’s “Bayesian Computation with R”. Note that while this might be homework, in my case it is not, as I am learning Bayesian methods in R because I think I might use it in my future analyses. Anyways, while this is a specific question, it probably involves basic understanding of Bayesian methods. So, in exercise 2.2, Jim Albert asks us to analyse the experiment of a penny throw. See here. We are to use a histogram prior, that is, divide the space of possible `p` values in 10 intervals of length `.1` and assign a prior probability to these. Since I know that the true probability will be `.5`, and I think it is highly unlikely that the universe has changed laws of probability or the penny is rugged, my priors are:               prior <- c(1,5,20,100,5000,5000,100,20,5,1)     prior <- prior/sum(prior)      along the interval midpoints               midpt <- seq(0.05, 0.95, by=0.1)      So far so good. Next, we spin the penny 20 times and record the number of successes (heads) and failures (tail). Easily done:               y <- rbinom(n=20,p=.5,size=1)     s <- sum(y==1)     f <- sum(y==0)      In my expeniment, `s == 7` and `f == 13`. Next comes the part that I do not understand: > Simulate from the posterior distribution by (1) computing the posterior > density of p on a grid of values on (0,1) and (2) taking a simulated sample > with replacement from the grid. (The function `histprior` and `sample` are > helpful in this computation). How have the interval probabilities changed on > the basis of your data? This is how it is done:               p <- seq(0,1, length=500)     post <- histprior(p,midpt,prior) * dbeta(p,s+1,f+1)     post <- post/sum(post)          ps <- sample(p, replace=TRUE, prob = post)      But **why do we do that**? We can easily obtain the posterior density by multiplying the prior with the appropriate likelihood, as done in line two of the block above. This is a plot of the posterior distribution: ![Plot of posterior distribution](http://i.stack.imgur.com/Kmz3V.jpg) As the posterior distribution is ordered, we can obtain results for the intervals introduced in the histogram prior by summarizing over elements of the posterior density:               post.vector <- vector()     post.vector[1] <- sum(post[p < 0.1])     post.vector[2] <- sum(post[p > 0.1 & p <= 0.2])     post.vector[3] <- sum(post[p > 0.2 & p <= 0.3])     post.vector[4] <- sum(post[p > 0.3 & p <= 0.4])     post.vector[5] <- sum(post[p > 0.4 & p <= 0.5])     post.vector[6] <- sum(post[p > 0.5 & p <= 0.6])     post.vector[7] <- sum(post[p > 0.6 & p <= 0.7])     post.vector[8] <- sum(post[p > 0.7 & p <= 0.8])     post.vector[9] <- sum(post[p > 0.8 & p <= 0.9])     post.vector[10] <- sum(post[p > 0.9 & p <= 1])      _(R experts might find a better way to create that vector. I guess it might have something to do with`sweep`?)_               round(cbind(midpt,prior,post.vector),3)                midpt prior post.vector      [1,]  0.05 0.000       0.000      [2,]  0.15 0.000       0.000      [3,]  0.25 0.002       0.003      [4,]  0.35 0.010       0.022      [5,]  0.45 0.488       0.737      [6,]  0.55 0.488       0.238      [7,]  0.65 0.010       0.001      [8,]  0.75 0.002       0.000      [9,]  0.85 0.000       0.000     [10,]  0.95 0.000       0.000      Furthermore, we have 500 draws from the posterior distribution that tell us nothing different. Here is a plot of the density of the simulated draws: ![enter image description here](http://i.stack.imgur.com/tGYqr.jpg) Now we use the simulated data to obtain probabilities for our intervals by counting what proportion of simulations are within the interval:               sim.vector <- vector()     sim.vector[1] <- length(ps[ps < 0.1])/length(ps)     sim.vector[2] <- length(ps[ps > 0.1 & ps <= 0.2])/length(ps)     sim.vector[3] <- length(ps[ps > 0.2 & ps <= 0.3])/length(ps)     sim.vector[4] <- length(ps[ps > 0.3 & ps <= 0.4])/length(ps)     sim.vector[5] <- length(ps[ps > 0.4 & ps <= 0.5])/length(ps)     sim.vector[6] <- length(ps[ps > 0.5 & ps <= 0.6])/length(ps)     sim.vector[7] <- length(ps[ps > 0.6 & ps <= 0.7])/length(ps)     sim.vector[8] <- length(ps[ps > 0.7 & ps <= 0.8])/length(ps)     sim.vector[9] <- length(ps[ps > 0.8 & ps <= 0.9])/length(ps)     sim.vector[10] <- length(ps[ps > 0.9 & ps <= 1])/length(ps)      _(Again: Is there a more efficient way to do this?)_ Summarize results:               round(cbind(midpt,prior,post.vector,sim.vector),3)                midpt prior post.vector sim.vector      [1,]  0.05 0.000       0.000      0.000      [2,]  0.15 0.000       0.000      0.000      [3,]  0.25 0.002       0.003      0.000      [4,]  0.35 0.010       0.022      0.026      [5,]  0.45 0.488       0.737      0.738      [6,]  0.55 0.488       0.238      0.236      [7,]  0.65 0.010       0.001      0.000      [8,]  0.75 0.002       0.000      0.000      [9,]  0.85 0.000       0.000      0.000     [10,]  0.95 0.000       0.000      0.000      It comes to no surprise that the simultion produces no other results than the posterior, on which it was based. Thus, **why did we draw those simulations in the first place**?


title: Panel-data exploratory data analysis

text: I have a dataset for around 40k firms over fiscal years 1950-2011 with about 430k firm-years. If I'm not mistaken I have panel-data. In addition, the firms are nested within 9 industries. I created a unique identifier `ticn` for each firm. Years are indicated by `fyear`. For now my variables of interest are yearly sales `sale`, yearly advertising `xad`, and yearly R&D expenses `xrd`. I have industry dummies indicated by `sicagg`. I am interested in the relationship between yearly sales and advertising/R&D expenditures, including some control variables. Currently I am in the exploratory phase of my research. _So my objective, I want to get a feel for the data, give descriptives and maybe make a few plots._ First I computed between and within descriptive statistics (mean, stdev, min and max). I also made scatter plots between sales, R&D and advertising. In addition, I plotted the time-series of the yearly average advertising expenses for each industry in a nice graph. Can you give me ideas for additional analyses? Thanks in advance!


title: T-test and the number of observations

text: I have a sample of the number of bank loans issued each year from the period 2004-2012. I have split the sample in two (2004-2008 and 2009-2012) and want to test if the number of bank loans issued in the second period (2009-2012) is significantly lower than the first period. Question:   When using the T-test, would N be the total number of bank loans (640 in this case) or the total number of years observed (9)? Thank you.


title: Can we find bounds on R-squared?

text: We know that as the number of independent variables increases, the coefficient of determination $R^2$ will increase but the adjusted $R^2$ may or may not increase. In the following question for the sake of simplicity I shall write only $R^2$ but it must be understood that the question applies to both $R^2$ and adjusted $R^2$. Further we shall make life easier by assuming all the conditions and assumptions of multiple regression are satisfied. **Question:** Consider a multiple regression where the the dependent variable $y$ depends on at most $n$ independent variables $x_1, x_2, \ldots, x_n$. For a given $k$, $1 \le k \le n$, we find best $k$ variable liner fit for $y$. Let us denote the coefficient of determination of this best fit by $R_{max}^2(k)$. Similarly we find the worst possible $k$ variable fit and we denote its coefficient of determination by $R_{min}^2(k)$. Trivially we have the following bounds $$ R_{min}^2(k) \ge R_{min}^2(1) $$ and $$ R_{max}^2(k) \le R_{max}^2(n) = R^2(n). $$ My question is can we find improve and express the above bounds in terms of non trivial functions involving $n$, $k$, $R_{min}^2(1)$ and $R_{max}^2(n)$. Is any additional assumptions is required to obtain such non trivial bounds? **Motivation** : I am currently working on linear modeling where I have a large number of independent variables and I need a way to determine how small or large the coefficient of determination will be for a given $k$. Currently I am following a various algorithmic approaches and writing a programs that gives above bounds. This method is not much useful because despite using the best known algorithms such as leaps computation takes a lot of time as the number of variables increases. Therefore I want to see if a theoretical bound is possible. **My progress so far** : Based on heuristic data I have generated using a computer program, I find that $R_{max}^2(k)$ approximately follows a logistic model $$ R_{max}^2(k) \approx \frac{R^2(n)}{1+ae^{-bk}} $$ where $a$ and $b$ are some local constants which depends on the data being analyzed.