title: Assessing stratum variability in surveys

text: What is the best approach to determining variability within different stratums in a survey ? The easiest approach would be to examine variability over any one response variable, but are there methods to take into account responses across a number of variables. The objective is to be able to adjust sample sizes for stratums with greater or lesser variability. Also all responses are categorical.


title: Variance structure with multiple covariates in GLS

text: I am building a GLS model following protocol in "Zuur, 2009. Mixed effects models..." on p.90. I have 5 continuous predictors. VarConstPower variance structure works best for me. At first the fixed part of the model includes all covariates. At that point I get the lowest AIC when variance part also includes all 5 covariates. But with model selection my fixed part shrinks to only 2 covariates. Can I still retain all 5 covariates in my variance part?


title: Analyzing a fractional factorial design

text: I designed an experiment using the fractional factorial design feature in Minitab. Can the results of the experiment be analyzed using something like linear regression and backward selection or is it best to use software that is aware that the experiment is a fractional factorial design experiment?


title: Help in understanding how to apply correctly False Discovery Rate adjustement

text: I am trying to understand how to correctly apply False Discovery Rate when comparing multiple hypothesis tests. Although I use here R code, my doubts are about the procedure, rather than the programming. I built a toy model in R made of 10000 hypothesis (e.g. genes expression) made over two 5 samples populations:               set.seed(620)     x = matrix(rnorm(10000*5),nrow=10000)     y = matrix(rnorm(10000*5),nrow=10000)      With these datasets `x` and `y` I know that all the null hypothesis are true. I now evaluate the p-values:               p = sapply(1:10000, function(i) t.test(x[i,],y[i,])$p.val)      As expected the number of p-values below 0.05 (or any other number) is 453, i.e. about 5% false positives as expected. Next I adjust the p-values using False Discovery Rate adjustment and estimate the q-values:               q = p.adjust(x, method = "fdr")      Now, if I understood correctly, selecting the hypothesis that have a q value of 0.05 one should get 5% of false discoveries (number of false positives divided by the number of discoveries). The number of hypothesis with q < 0.05 is 0. I think that it might be because, since all the null hypothesis are true, no matter how I choose q the false discoveries will always be 100% among the discoveries (this is also how I explain to myself that most of the q-values are close to 1). Next I substituted the last one hundred rows of y with numbers sampled by a normal distribution with mean 3, and estimate the p and the q-values:               y[9901:10000,] = rnorm(500, mean=3)      p = sapply(1:10000, function(i) t.test(x[i,],y[i,])$p.val)     q = p.adjust(p, method = "fdr")      After these modifications, the number of p-values < 0.05 increases to 544 and 98 out of the 100 hypothesis that should be rejected are detected. However the number of hypothesis having q-values < 0.05 is surprisingly low: only 9. They are all hypothesis that should be rejected, so it seems to me that the false discovery rate has been kept to 0 rather than to 0.05. If I accept, for example, the hypothesis with q-value = 0.5 I end up with accepting 95 hypothesis. Of these 95, 67 are true discoveries and 28 are false discoveries. Therefore the FDR is of 28/95 = 0.3 and not 0.5 as I would expect. Is there something I have not understood properly? Why are the result I get so different from the one I would expect theoretically?


title: One-way repeated measures ANOVA with G*Power

text: I have made an experiment with 28 subjects. All made the experiments under 3 different conditions (A, B and C). So I have a table like this:               No. A B C     1   3 6 4     2   6 2 6     ...      Now I just want to perform a RM-ANOVA or a Friedman test or a Quade test. That's not the problem. I want to use G*Power to determine the necessary alpha value to detect e. g. a "medium" effect. I do not get how to use G*Power in this case. Here is a screenshot of how I do it for a t-test: ![enter image description here](http://i.stack.imgur.com/6mWCJ.png) Now I want to do the same for more than 2 conditions. Can somebody point out how to do that? Here's a screenshot of the possibilities: ![enter image description here](http://i.stack.imgur.com/GmjhZ.png) **EDIT:** I looked it up, and apparently it is called _One-way repeated measures ANOVA_. So how can I to _One-way repeated measures ANOVA_ with G*Power to determine the required alpha value (given effect size and so son)?


title: Standard deviation of the sum of two normally distributed random variables

text: $X\sim N(52,6)$, $Y\sim (40,8)$. What's the standard deviation of $Z=X+Y$? I'm considering to transform the linear relationship to matrix form $$Z=\begin{pmatrix} 1& 1\\\ \end{pmatrix}\begin{pmatrix} X\\\ Y \end{pmatrix}$$ and apply $$\Sigma_{z}=A \Sigma^{-1} A'$$ However, I don't know the next step. How to calculate $\Sigma$? Thank you!


title: Application of Bayes' theorem

text: **Background:** My professor gave us several copies of his hand-writing handouts last week. Due to the low quality of printer, there are several words in the handouts cannot be recognized. However, I found that I can " _guess_ " what the professor has written on the handouts by reading context carefully. **Question:**   1. Why could we _guess_ and recognize vague words on the handouts?   2. Is it related to statistics?    3. If we find out how statistics works in the recognition process, can we apply it to other fields? **My idea:** _**Assumption_** :   * Define $0$ represents the recognizable degree we cannot recognize the words at all.   * Define $100$ represents the recognizable degree we can recognize the words easily. _**A simple example:_**                                        I have ki_ _ed her.      In this sentence, "I","have","ki ed"and"her" is easy to recognize and they're 100. "_ _" represent two letters that we can not recognize and they're 0.   1. If we read the sentence in a romantic novel, it is reasonable to _guess_ the two missing letters are "ss",that is, " I have kissed her."   2. If we read the sentence in a detective novel, it is reasonable to _guess_ the two missing letters are "ll",that is, "I have killed her." **Bayes' theorem:** $P(ll \mid D)$= the probability that we find "kill" in the context of detective novel. $P(ll) can be simply considered as the proportion of "kill" in a detective or a romantic novel. $$P(ll \mid D)=\frac{P(D \mid ll)P(ll)}{P(D)}$$ $$P(ss \mid D)=\frac{P(D \mid ll)P(ll)}{P(D)}$$ Here $P(D)$ is a constant. $$P(ll \mid D) \propto P(D \mid ll)P(ll)$$ That's it ! It's simple and natural. _**What's the relationship between the background and the example_** The recognizable degree of vague words range from 0 to 100. The example,two of the letters are 0 and others are 100, is a extreme case. **My question:**   1. Is there any mistakes or problem?   2. Could anyone provide any useful resources like articles, books etc to help me to go further on this topic?   3. I think is important to prove that human beings can naturally and automatically find $P(ll \mid D) \propto P(D \mid ll)P(ll)$ in their brains. How can I use statistics tools to prove that? Design an experiment? * * * _**Extension_** : Why could we read the words on the image below? How does our brain automatically correct the order of letters? I guess Bayes helps again. ![enter image description here](http://i.stack.imgur.com/MdT8Y.jpg) _**Extension2_** : For many Chinese students, they have never learn traditional Chinese. However they can read a newspaper written by traditional Chinese easily. When they are asked to recognize a single tranditional Chinese character, they failure in most cases. Why? **Thank you!**


title: Choosing the best analysis

text: I need to do a data analysis in R. My data consists of 2 populations of animals that have 2 different types of genomes. The children are always born as twins. I also know the mothers' weight and the combined weight of the children after 14 days. I want to fit a model that explains the variability of the twins' weights after 14 days. Data:               population 1 and 2     genome 1 and 2     twins both male, both female, one female and one male     mother's weight     twins weight after 14 days      I'm planning to do an ANCOVA, but is there a better analysis I should use? All suggestions are welcome. Big thanks for the help.


title: Can a neural network output represent a posterior probability?

text: I seem to remember from years ago when I first read Bishop's ANN book that it is possible to construct a neural network such that the outputs should represent the posterior probability that I would have found if I used the likelihood approach to separate the same data. Is this true or have I made this up in my mind? If it is true what are the conditions, I seem to remember it has to have no loops and limits on the number of hidden layers. Although I also seem to remember a result from Komogorov showing that many hidden layers aren't helpful if the network is well trained. Unfortunately I can't find the book anymore and I can't seem to find much on this result which makes me worried it's me. I think I probably need to flesh out my question a bit as pointed out by @bayerj I have a few variables (I'll assume I only have one for now) which I am using to find the probability of one of three hypotheses being true given the three variables. I have some data which which a can construct PDFs. I am using $$P(H_{k}|x) = \frac{P(x|H_{k})P(H_k)}{P({x})}$$ where P(x) is the unconditional probablity density \begin{equation} P(x) = \sum P(X|H_k)P(H_k) \end{equation} I'm trying to verify whether it is possible to train a network, using the data with which I construct my PDFs, such that the value of the output node equals P(H_k|x). I'm sure this is possible I just can't remember the working. Since I posted the question I have found http://www- vis.lbl.gov/~romano/mlgroup/papers/neural-networks-survey.pdf which I think is saying what I want to do is possible.


title: Distribution fitting on empirical PDF?

text: I have a class of distributions that I think may have generated some data I have, and would like to find the best fit distribution and corresponding parameters. But the data I have aren't samples from a distribution, they are an empirical PDF.


title: Support Vector Machines - Kernel Functions/Soft Margin SVM

text: I had these questions in an exam today. State True or False and explain.   1. If k1(.,.) and k2(.,.) are two valid kernel functions, then if h = k1 - k2, is h(.,.) a valid kernel function?   2. A standard soft margin SVM is used to classify data set. We have a fixed C parameter. Two different algorithms A1 and A2 are used to obtain the support vector set {$S: \alpha_i > 0$}. Call them S1 and S2. Is S1 = S2 in all cases? Assume both algorithm use the same kernel function. I had no clue about these questions. I only knew about kernel functions, but not about operations between kernel functions. I only guessed: as kernel function need to be positive semi definite (PSD), the difference between two kernel functions need not be PSD. Hence FALSE. For the second one, I guessed as $\alpha_i$ can be different among the two algorithms, the number of support vectors can differ as well. Hence FALSE again.


title: Correlation with non-normal distribution

text: I want to find the correlation between a continous variable with a non normal distribution, positively skewed, and a categorical variable with 4 categories. The only test I can find that fits these variable is the kruskal wallis test, but I am not sure how I can use this test to find the correlation or relationship between my two variables. Can anyone suggest how this relates to correlation or another test to use with my data? Rookie with statistics, so please let me know if any more information is required.


title: Advantage of using ANOVA over a kruskal wallis test

text: This question is a follow on from this question - Correlation with non-normal distribution. An ANOVA explanation was given here - Correlation between vegetation and erosion \- that could be applicable to my question. Further to this, I was hoping someone could tell me if there was an advantage to using ANOVA over a kruskal wallis test, and if ANOVA can be used in this way to find the correlation between a categorical varaible and a continuous variable with a non-normal distribution.


title: How to choose between two competing linear models?

text: I have the following data:               t       mean     147     1.4     143     3     137.5       1.8     133     1.9     129.5       1.8     124.5       2.5     115.5       1.9     107     2.5     102.5       6.3     98.5        6.5     94.5        5     89      5.5     81      4.8     73      9.3      To me, the slope looks more exponential than linear when plotted as a scatterplot. I've been using the following code in R:               data<-read.csv("regression.csv")     attach(data)     plot(t,mean)     data.lm<-lm(mean~t,data=data)     summary(data.lm)     data.exp<-lm(log(mean) ~ log(t) ,data=data)     summary (data.exp)     AIC(data.lm, k=2)     AIC(data.exp, k=2)      data.exp, the exponential regression, has a much lower p-value and a much lower AIC score than data.lm, the linear model: 6.869e-05 vs. 0.000194, and 11.4641 vs. 52.22926. But (how) can I demonstrate that the data fits an exponential line better than a straight line? Is the use of AIC legitimate here? Sorry to ask such as simple question but I've looked online and haven't found an answer. Thank you!


title: A robust R package to do MCMC and Gibbs sampling

text: I need to make linear model for which I need to do Gibbs sampling in MCMC simulations. The model needed to be fitted is a linear mixed model. Please suggest me for a robust R package for this task.


title: subtract beta coefficients two models

text: I have two regression models and I want to subtract one beta coefficient in model 1 from model2. Obviously I can just do this on a calculator, but does anyone know of a test I could use in R to test the difference for me, and give me back confidence intervals or standard errors. I've tried running a t.test:               t.test(betamod1,betamod2)        and get this error               not enough 'x' observations      I have also tried trying to calculate a z-score using               z = (betamod1-betamod2) / sqrt(stanerrbetamod1^2 + stanerrbetamod2^2)      Thank you in advance


title: Question on expected value of an estimator

text: I have been given the pdf: $$f_{\theta}(x)=\frac{x}{\theta^2}e^{-\frac{x^2}{2\theta^2}}$$ and I want to know if the MLE estimator $\theta$ is unbiased. Attempt: I want to maximize the function so I take the natural logarithm and then I differentiate with respect to $\theta$ as follows: $$\frac{d}{d\theta}ln(f_\theta (x))=\frac{-2}{\theta}+\frac{x^2}{\theta^3}$$ which equals zero iff $\theta =x^2/4$. Thus $\hat{\theta}=\frac{x^2}{4}$ Now I want to evaluate $E[\hat{\theta}]$ which I'm having trouble to implement. Is it $$E[\hat{\theta}]=E[x^2/4]=\int_0^{\infty}\frac{x^2}{4}f_\theta (x^2/4)dx$$ or am I misunderstanding the definition of $E[\hat{\theta}]$? Edit:differentiation fixed


title: Binary classification of DNA sequence motifs

text: [I apologize in advance for the length of this question. but it seems unavoidable.] I'm trying to do binary classification, using Bayesian methods, and I'm having some trouble figuring what methods I need to do so. Bottom line - I need some kind of decision rule. I've tentatively settled on a version of the Bayes factor, outlined below, but am open to other suggestions. I know little about Bayesian statistics or decision theory, and this project has become to some extent an exercise in learning some Bayesian theory. My formal statistical education has been rather frequentist in nature. I also know a little about machine learning, but that doesn't seem relevant here. Here is a description of the background. The data is a set of sequences of letters of equal length, say $n$. The letters are from some alphabet with $r$ letters. The specific application is for DNA sequences i.e. strings composed from the four letters A, C, G, T. So here, $r=4$. The aim is to use the elements of such a set of DNA sequences (called motifs) to predict which other sequences belong to the family. Personal/historical note: this was started some years ago as a project in motif finding, and was abandoned. I've recently picked it up again. The approach outlined here is outside the mainstream of motif finding, as far as I can tell, because motif finding is mostly concerned with finding new motif families, not extending existing ones. This approach could be used to find new members of existing families, but I don't see how it could be used to find new motif families. Let $\mathbf{X}=(X_1,\dots, X_n)$ be a random vector. We model these sequences as a set of i.i.d. random vectors, $\\{\mathbf{X^{(i)}}=(X_1^{(i)},\dots, X_n^{(i)})$, $i\in \mathcal{I}\\}$, where $\mathbf{X^{(i)}}\sim\mathbf{X}$ $\forall i\in \mathcal{I}$. We define a family of probability distributions for X. We first try to find the member of this family which is the best fit for the data, and then use it for prediction. Our assumption is that $\mathbf{X}$ can be divided into disjoint subsets, such that variables in each subset are dependent among themselves, and the subsets are mutually independent. It suffices to define the joint probability distribution on one of these subsets. Let $\mathbf{S}=(X_{j_1},\dots, X_{j_t})$ be a sub-vector of $\mathbf{X}$ for some subset $\\{j_1,\dots, j_t\\} \in \\{1,2\dots, n\\}$. So, $\mathbf{S}$ takes $k=r^t$ possible values. Then the joint probability distribution is defined as a categorical probability distribution on $\mathbf{S}$. \begin{align}\label{prob} P(S=(s_{1}, s_{2}, \dots, s_{t})) = p_{s_{1}s_{2}\dots s_{t}} \end{align} So the distribution of $\mathbf{S}$ is defined by $k$ probabilities. Observe that the random variables $\\{X_{j_1},\dots, X_{j_t}\\}$ are not independent for general values of this probabilities. Let us restrict our attention to the sequence of subsets of $\\{\mathbf{X^{(i)}}, i\in \mathcal{I}\\}$, corresponding to $\mathbf{S}$, which we call $\mathbf{Y} = \\{\mathbf{S}^{(i)}, i\in \mathcal{I}\\}$. Observe that this is a sequence of i.i.d. categorical variables. Then the likelihood of $\mathbf{Y}$ is \begin{align*} P(\mathbf{Y}|\mathbf{p}, \mathcal{M}) = \prod_{i=1}^k p_i^{f_i} \end{align*} where $(f_1, \dots, f_k)$ are the frequencies of the $(s_{1}, s_{2}, \dots, s_{t})$ in Equation~\ref{prob} for some fixed ordering, and $\mathcal{M}$ is the chosen model, namely the division of $\mathbf{X}$ into subsets. The corresponding uninformative prior for $\mathbf{p}$ is \begin{align*} \pi(\mathbf{p}) = (k-1)! I(\textstyle\sum_i p_i = 1) \end{align*} This is a special case of the Dirichlet distribution \begin{align*} \pi(\mathbf{p}) = \frac{\Gamma(\sum_i{\alpha_i})}{\prod_i\Gamma(\alpha_i)} \prod_{i\in\mathcal{I}} p_i^{\alpha_i -1} I ( \textstyle\sum_i p_i = 1 ) \end{align*} with $\alpha_i=1$ $\forall i$. Now, one can use Bayesian model selection to select a suitable model, by searching for the model that maximize $P(\mathbf{Y}|\mathcal{M})$. Let us assume we have found a suitable candidate for this model, namely $\mathcal{M}_1$. Call the corresponding uncorrelated model $\mathcal{M}_0$. The idea is to try to discriminate between $\mathcal{H}_1=\mathcal{M}_1$ and $\mathcal{H}_0=\mathcal{M}_0$, to try to use this model for prediction. To test this, we use k-fold cross-validation. We divide a data set of related sequences (motifs), into $K$ equal parts. We use $K-1$ of these parts as a training set, and the remaining part as a testing set. Then we apply a decision rule to every element in the testing set to determine whether the element is a better fit for $\mathcal{M}_1$ or $\mathcal{M}_0$. The obvious approach is to use the Bayes factor, namely \begin{align}\label{bayesfact} K = \frac{P(D|\mathcal{M}_1)} {P(D|\mathcal{M}_0)} = \frac{\int P(D|\theta, \mathcal{M}_1)P(\theta|\mathcal{M}_1)\, d\theta} {\int P(D|\theta, \mathcal{M}_0)P(\theta|\mathcal{M}_0)\, d\theta} \end{align} where $D$ is an individual element in the testing set. However, it turns out that with an uninformative prior, the values of this are the same for all sequences, since, with an uninformative prior, there is not enough information to distinguish between the sequences. After some confusion, it became clear that one needs an informative prior here, otherwise one cannot do anything. It then occurred to me that one can use the posterior distribution for the parameters instead of the prior, namely $P(\theta|\mathcal{M}_1, D_0)$ in the place of $P(\theta|\mathcal{M}_1)$, in Equation~\ref{bayesfact}. Here $D_0$ corresponds to some other data, that is used to get an informative distribution for the parameters. In this case, the obvious choice is the training set. I wasn't sure if that was actually done or not; however after some searching, I came across this variation described as pseudo-Bayes factors in Ch 5 of Robert's ``Bayesian Choice''. However, from his description (I have not looked into this otherwise) the main motivation for doing this is that bane of Bayesians, improper priors. I.e. if one conditions on a training sample. the posterior corresponding to an improper prior may be proper. I don't know if anyone is using this in the context I have described. So, let $K$ be a pseudo-Bayes factor. Then \begin{align*} K = \frac{P(D|\mathcal{M}_1, D_0)} {P(D|\mathcal{M}_0, D_0)} = \frac{\int P(D|\theta, \mathcal{M}_1, D_0)P(\theta|\mathcal{M}_1, D_0)\, d\theta} {\int P(D|\theta, \mathcal{M}_0, D_0)P(\theta|\mathcal{M}_0, D_0)\, d\theta} \end{align*} where $D_0$ is the training set, and $D$ is an individual element in the testing set as before. This is the same as the ratio of the the post-predictive distributions corresponding to $\mathcal{M}_1$ and $\mathcal{M}_0$ respectively. As usual for Bayes factors, greater than 1 accepts the alternative hypothesis, less than 1 accepts the null. What I did for testing was permute the testing set, which is an array, along the columns. This destroyed any correlation structure that might be present. Call $D^{\prime}$ the testing set, and $D^{\prime}_p$ the permuted set. Then I took all the sequences from $D^{\prime}$ and $D^{\prime}_p$ (which are two sets of equal size) and calculated the log of the pseudo-Bayes factor for them. One would expect the sequences in $D^{\prime}_p$ to fall to the left of 0 and $D^{\prime}$ to the right of 0. Preliminary testing gave something close to that. I was wondering if this was a reasonable testing framework, and how to summarize the results. The motivation for the permutation was to get a dataset which had a similar (empirical) marginal distribution to the real dataset, but also which had had any possible correlation structure removed from it. The idea was to construct a set that would be as hard to distinguish from the original dataset as possible. It would nice to have a more formal criteria than this, though. So, here are some questions.   1. Is it reasonable to test on $D^{\prime}$ and $D^{\prime}_p$? Can one give a formal reason for doing this?   2. I thought I would plot the PBF for one instance of $D^{\prime}$ and $D^{\prime}_p$. What would be a good way to graphically represent this? Some kind of bar plot?   3. I thought of summarizing the overall classification results by using a classification table, i.e. a 4-way table with the counts that     * correctly accept $\mathcal{H}_0$     * correctly accept $\mathcal{H}_1$     * incorrectly accept $\mathcal{H}_0$     * incorrectly accept $\mathcal{H}_1$ Does this look reasonable? It is Ok to add total classification counts across all training sets? That seems like a good summary approach.   4. Does this approach using a pseudo-Bayes factor makes sense, and are there any good alternatives to doing it this way?   5. One is expected to show how this method performs in comparison with other motif finding approaches. Do you have specific suggestions about which methods to compare it to?   6. There is something called the 0.632+ Bootstrap, which is supposed to be an improvement on regular cross-validation. Is this something worth considering? I'm also a little fuzzy about the justification of the Bayes factor, but that can wait for another question. If you think the notation in this question could be improved, please comment. Thanks.


title: Bayesian analysis problem

text: The problem in hand is that the prior distribution which I have received from experts (loan recovery data) ranges from 0 to 100%. Thus a beta distribution was assumed. Where as the actual data shows that loan recovery can be more than 100% due to fees and interest charges. Thus the likelihood function does not follow a beta distribution simply because the values may be greater than 100%. Trying to use conjugate prior method to combine distributions to arrive at the posterior distribution. First I thought the beta-beta conjugate prior method will suffice. Rescaling the likelihood function to 0% - 100% does not seem correct. What do I do. Please help


title: Do null and alternative hypotheses have be to exhaustive or not?

text: I saw a lot of times claims that they have to be exhaustive (the examples in such books were always set in such way, that they were indeed), on the other hand I also saw a lot of times books stating they should be exclusive (for example $\mathrm{H}_{0}$ as $\mu_1=\mu_2$ and $\mathrm{H}_{1}$ as $\mu_1>\mu_2$) without clarifying the exhaustive issue. Only before typing in this question I found somewhat stronger statement on the Wikipedia page \-- "The alternative need not be the logical negation of the null hypothesis". Could someone more experienced explain which is true, and I would be grateful for shedding some light on the (historical?) reasons for such difference (the books were written by statisticians after all, i.e. scientists, not philosophers). Thank you in advance.


title: Choosing variables for Discriminant Analysis

text: I've 110 variables & 200 data points. Of this 110 variables, one is group variable (say "brown eye","blue eye"). I want to use discriminant analysis to classify the groups based on remaining 119 variables. Since the variables are large, to get a meaningful result I need to reduce the number of variables. So, the are 3 options to me: 1) Stepwise Discriminant Analysis: I don't want to use this method as I'm biased towards it. 2) Classification Tree Method: This method will give an idea about which variables affecting the eye color. Since the dataset is small, I'm apprehensive of using this method. 3) Principal Component Method: This method I can use. But I prefer to keep the original variables. My question is can anybody please suggest me some other method to select variables for discriminant analysis.


title: How to show that the variances of 2 sets of 3D points are different?

text: I have 2 point clouds (3D points). I can visually tell that the spread in one cloud is much larger than the other, and I have also plotted their error ellipsoids. Now, I'm looking for a statistical test which proves that the the spread is larger in one cloud. So my question is, if I use a statistical test like Levene/Brown-Forsyth/Box-M to test if the variances of two populations are not equal, does it mean that I have statistically shown that the variance (or the spread) or the clouds are different?


title: How to compute theoretical compression limit?

text: Assume we have a sensor field with dimension M*M. In order to apply any data compression technique, first I want to know what is the compression limit or minimum entropy of the entire sensor field. How could I compute the minimum entropy or compression limit for the sensor field? or Actually I want to have the theoretical compression limit. Let's put the problem for an image. I want to know whether there are any mathematical methods to calculate the theoretical compression limit. Please let me know or suggest any readings to formulate the problem. Thanks


title: Varying group coefficients in lme4

text: All, I am estimating a multilevel logistic regression with group predictors, but am unclear about some of the advice given by Gelman and Hill (2007) in their book. Therein, they recommend allowing every coefficient to possibly vary, given a large enough N. Does that include group predictors as well? They weren't clear, treating "varying slope" as just another complexity you can incorporate into a mixed effects model in lme4 along with group predictors (see: p. 549 in their book). For example, I have roughly 50,000 observations with a binary response (plenty large N). Predictors exist at two levels, such that my model looks like: `M1 <- lmer(Y ~ X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 + (1 | group), family=binomial(link="logit"))` X1:X4 are individual-level predictors and G1:G4 are group-level predictors, thus: a multilevel model. Does their recommendation of treating all coefficients as potentially variable mean including even the group predictors within the random effect, such that: `M2 <- lmer(Y ~ X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 + (1 + X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 | group), family=binomial(link="logit"))` I ran M2 and it gave sensible estimates. AIC/BIC suggest much better fit than M1. I'm just unsure if it's appropriate since, unlike individual-level predictors, the group-level predictors are not going to vary in a given group. It will obviously vary across groups, though. Further, if this is not an incorrect way to approach it, how suspicious should I be if one of the group predictors of interest is statistically insignificant as a stand-alone fixed effect (varying intercept model like M1), but is significant as a fixed effect in a varying slope model like M2? Thanks for any input and feedback on this topic. I really appreciate it.


title: What machine learning techniques are especially susceptible to "over-tuning" of their hyperparameters?

text: For example, the random forest algorithm is not especially susceptible to "over-tuning" because it has 1 hyperparameter, "mtry," and mtry typically does not have a large influence on the result of the algorithm. On the other hand, I imagine the "knn" algorithm is prone to "over-tuning", as it's hyperparameter "number of neighbors" has a wide range of possible values, and exerts great influence on the result of the algorithm.


title: How to select the exponential decay constant for weighting in proc logistic?

text: I am trying to predict the sales conversion using proc logistics in SAS. Right now I have around 3 months of data, and it will gradually grow to more than an year over time. My intuition is that the most recent data is a better predictor of the event, and hence I want to weight the data using the exponential weighting method. For this I have two questions:   1. How do I select an optimum decay constant for the exponential weighting? Is there a scientific way to go about it instead of just trial and error?   2. As my data grows larger, is there a way to ensure that exponential weighting changes or remains optimized? This may also cause overfitting. I tried searching this on the forum but could not find an explanation in this context. Any help will be much appreciated. Thanks.


title: Exposure in Negative Binomial regression and other distributions

text: For a Poisson regression we can model exposure $\epsilon_i$ in observation $Y_i$ as $Y_i \sim Poisson(\epsilon_i*\lambda)$. For instance, in a Poisson regression, if we observe: $y = \begin{bmatrix} 2 & 3 & 4 & 5 \end{bmatrix}$ With exposures $\epsilon = \begin{bmatrix} 2 & 1 & 1 & 1 \end{bmatrix}$ And covariate $x = \begin{bmatrix} 2.1 & 3.1 & 4.3 & 5.2 \end{bmatrix}$ We'll have a regression with the same coeficients as if we have observed: $y = \begin{bmatrix} 1 & 1 & 3 & 4 & 5 \end{bmatrix}$ With exposures $\epsilon = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 \end{bmatrix}$ And covariate $x = \begin{bmatrix} 2.1 & 2.1 & 3.1 & 4.3 & 5.2 \end{bmatrix}$ This because is Poisson can be written as (proportional to) $\lambda^{(\sum Y)} e^{(\sum \epsilon)}$, so it doesn't matter as long as we keep the total sum of y's and exposure constant and this, correct me if I'm wrong, is a tautology of the property of the Poisson mean and variance being equal. This does not follow for a negative binomial because its variance is different from its mean. Here's an example in R to illustrate that this follows for Poisson, but it doesn't for negative binomial:               require(MASS)     #These 3 have the same coefficients:     glm(c(Days[1]*2,Days[-1]) ~ Sex/(Age + Eth * Lrn) + offset(log(c(2,rep(1,nrow(quine)-1)))), data = quine, family=poisson())     glm(Days ~ Sex/(Age + Eth * Lrn), data = quine, weight=c(2,rep(1,nrow(quine)-1)), family=poisson())     glm(Days ~ Sex/(Age + Eth * Lrn), data = rbind(quine[1,],quine), family=poisson())               glm.nb(c(Days[1]*2,Days[-1]) ~ Sex/(Age + Eth * Lrn) + offset(log(c(2,rep(1,nrow(quine)-1)))), data = quine)     #the above produces different coefficients from the following two:     glm.nb(Days ~ Sex/(Age + Eth * Lrn), data = quine, weight=c(2,rep(1,nrow(quine)-1)))     glm.nb(Days ~ Sex/(Age + Eth * Lrn), data = rbind(quine[1,],quine))      So, can we still use exposure on it in the same fashion like: $Y_i \sim Negative\\_binomial(\epsilon_i*\mu, \phi)$ (parametrized using mean and overdispersion)? What about other distributions, does this idea of exposure could be applied, say, for a binomial? How? We would be obligated to parametrize it using the mean?


title: Test for measuring change when you have same primary sampling units but different respondents in two waves of longitudinal survey

text: Which test should I apply for measuring change in proportions and means when I have same primary sampling units but different respondents in two waves of longitudinal survey? Should I apply McNemar's test for proportion and dependent-samples t test for means, or which test to study the difference between two waves? How can I do this in STATA?


title: Probability to draw a black ball in a set of black and white balls with mixed replacement conditions

text: When a black ball is drawn, it is not replaced in the set, but white balls are replaced. I have thought of this, with the notations:   * $b$, $w$ the initial number of black and white balls   * $x_i = (b - i)/(b + w - i)$  The probability of drawing a black ball $Pb(n)$ after n draws: $$\eqalign{ Pb(0) &= x_0\\\ Pb(1) &= (1-x_0)x_0 + x_0x_1\\\ Pb(2) &= (1-x_0)^2x_0 + x_0x_1(1-x_0)+ x_0x_1(1-x_1) + x_0x_1x_2 \\\ Pb(n) &= \sum\limits_{k=0}^{n-1} (\prod\limits_{i=0}^k x_i \prod\limits_{i<=k}^{n-k\ terms} 1-x_i) }$$ This sum seems infinite with n, even if some terms are null since $x_{i \ge b}=0$ Except $b=1$:   $Pb(n) = (1-x_0)^nx_0 $ For $b=2$:   $Pb(n)= x_0(1-x_1)^n + x_0x_1\sum\limits_{i+j=n-1} (1-x_0)^i(1-x_1)^j$ Is there a known solution to this problem?


title: How to analyze unequal samples in linear and multiple regression

text: I have three variables (sample size is mentioned below), and I want to analyze data using regression models as recommended by Judd and Kenny (1981) to see if $cc$ mediates the relationship between $spdm$ and $cty$ (I'm using abbreviations; see below for variable names). I have a unequal sample size problem, and I don't know how to handle it. I have read a few books and went through many articles, but none have mentioned or analyzed samples of unequal sizes. **I want to run 3 models, of which 2 are with unequal samples:** model 1: $spdm = \beta \cdot cty + e$ model 2: $cty = \beta_1 \cdot cc + \beta_2 \cdot spdm + e$ (where $e$ is error) **My Questions:**   * How can I analyze my data in regression? Is it possible to analyze samples of unequal size?   * I read somewhere on the web that you can draw a random sample from the big sample to make it equal to other samples. Can I do that? _Please provide the source of any book or article._ **Sample sizes and variables:**   * a = 100 ($spdm =$ supervisor's participation in decision making)   * b = 100 ($cc =$ creative climate)   * c = 500 ($cty =$ creativity)


title: Manually fitting a mixture distribution in matlab

text: I am trying to fit a mixture model containing a gamma and an exponential distribution: The general form, using the pdfs, is: p * gammapdf + (1-p) * exponentialpdf. The pdfs for the Gamma and Exponential, respectively are as follows: ![Gamma PDF](http://i.stack.imgur.com/iVv1O.png) ![Exponential PDF](http://i.stack.imgur.com/WBIZL.png) In matlab I have coded the mixture of these two as follows:               (p * ((1/(gamma(k) * theta^k)) * (x(i)^(k-1)) * (exp(-x(i)/theta)))) + ((1-p) * (1/lambda) * exp(-x(i)/lambda))      where p is the probability of the gamma and 1-p is the probability of the exponential, gamma is the internal matlab call to the gamma function, k is the shape parameter, theta the scale parameter of the gamma and lambda the rate of the exponential. x is the vector containing all of the data and i indexes x (I'm using a for loop to evaluate each data point with the model, I could vectorize it, but the speed is actually pretty good, all things considered). However, the predicted data are no where close to what one would expect. My intuition is that there is something wrong with how I coded the gamma portion of the mixture distribution. Upon looking at this, does anyone have any suggestions for how I might improve this? Additionally, if you have time, a second mixture distribution I am trying to manually fit is a mixture of the exponential and inverse gaussian distributions. Below is the pdf for the inverse gaussian ![Inverse Gaussian](http://i.stack.imgur.com/gdbwr.png) And the code for this mixture model:               (p * ((lambda/(2*pi*(x(i)^3)))^0.5 * (exp((-lambda*((x(i) - mu)^2))/(2*(mu^2)*x(i)))))) + ((1-p) * (1/lambda_exp) * exp(-x(i)/lambda_exp))      Where p is the same as above, except p corresponds to the inverse gaussian, lambda is the mean of the inverse gaussian, mu the standard deviation, and lambda_exp is the rate of the exponential. And again, x is as above. I am using maximum likelihood estimate, where I take the natural log of each value and then sum those values to get a log-likelihood (LL) value that is then fed in simulated annealing (SA) minimization function in matlab to find the best parameter values. Note, that since SA minimizes, I put a negative sign in front of the LL value, which effectively causes SA to maximize. I broke out the individual components of my mixtures and have confirmed that using these process I get the same parameter values when fitting the single distribution as the built in functions in matlab. I'm now thinking the issue is in how I get the predicted values and then plot them. My tactic here was to take the same equations from above and feed in the parameter values and the data to get predicted values. This however and obviously, gives me the pdf of the data and not actual predicted values. My question then is, given that I have a mixture of two distribution pdfs, how do I get actual predicted values from them so I can plot them against each other?Thanks.


title: Why is 30 touted as a sample size when the tables say otherwise?

text: If I have a population of 5000, statistics suggests that I need to sample about 350+ to get a confidence interval of 95% with margin of error 5%. So why do I see that sometimes we can get away with as low as 30?! That seems so much lower than what all these tables and calculators suggest.


title: Understanding coxph output in R

text: I am attempting to fit a Cox proportional hazard model to my data. I think I have the formula correct but am having trouble understanding the output. I have tried looked through the documentation and it is hard for me to understand. Any help would be greatly appreciated, thank you! The formula:               coxfit1 <- coxph(Surv(days, status)~GENE1, data=dataset1)     summary(coxfit1)      Where "days" is days until an event occurred (or last known followup if no event), "status" is an event (recurrence), GENE1 is expression data of a gene that I am testing if it has an effect on recurrence. The output:               Call:     coxph(formula = Surv(days, status) ~ GENE1, data = dataset1)          n= 34, number of events= 22                       coef exp(coef) se(coef)     z Pr(>|z|)        GENE1    0.6370    1.8908   0.2362 2.697  0.00699 **     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                    exp(coef) exp(-coef) lower .95 upper .95     GENE1        1.891     0.5289      1.19     3.004          Concordance= 0.618  (se = 0.068 )     Rsquare= 0.166   (max possible= 0.98 )     Likelihood ratio test= 6.17  on 1 df,   p=0.01298     Wald test            = 7.27  on 1 df,   p=0.006993     Score (logrank) test = 7.81  on 1 df,   p=0.005198      Now, this is one that is obviously significant, but what do the different parts of this output mean? Where is the hazard ratio??? And which of this information is appropriate for reporting?


title: Estimating confidence of a prediction

text: Given a set of features vectors $X=\\{\vec{x}_1,..,\vec{x}_n\\}$, binary ground truth data $Y=\\{y_1,..,y_n\\}$ and continuous prediction $\bar{Y} = \\{\bar{y}_1,..,\bar{y}_n\\}\in [0,1]$, I want to perform some regression analysis to estimate a function $f(\vec{x}_p) \approx e_p$ where $e_p = abs(y_p - \bar{y_p})$. In practice I want to be able to predict the confidence of a classification based on some external features that do not belong to the classifier itself (that is treated as a black box). An example could be: I have a classifier that predicts whether a fish is a salmon or a carp. Given a set of features such as weight,size,eye color etc, I want to find a measure that tells me how reliable the classifier is based on those features. E.g. I might find out the classifier is unreliable for fishes of size smaller than 30cm. I am not an expert in the field and I am looking for the best way to proceed.


title: Calculus resources with a probability and statistics slant?

text: I'm taking a statistical learning class, and I need to brush up on my calculus. I already took calculus and probability several years ago, but I have a bad memory. Are there any good books or general resources that give a review of calculus with a probability and statistics slant? I noticed that most calculus books give examples in the context of physics. In general, I'm looking for a book like Calculus Refresher, except more targeted towards probability and statistics.


title: How to validate a Multinomial Logit and Probit Model fit?

text: I would like to know how do you determine the performance of your models. That is, if you fit a multinomial logit or probit model for un-ordered discrete choice. What do you use to evaluate whether you have a good model? Please provide me with any reference material that I could read more on this topic. Websites, articles or book references on this topic will be a great help. I am kind of stuck! I have R, SAS, Stata, SPSS, and Minitab available.


title: What to do if samples aren't comparable?

text: I have an unbalanced panel data set, which covers a set of about 10 variables I am very interested in as controls in a regression. However, only a few of these variables are available in all waves: most of them are only available for a few waves down to two which are only available for 5 waves. Additionally, these variables have high rates of missingness. If I just create a sample where all variables are available I lose a massive share of observations and years. So I split the sample in several subsamples each depending on the availability of these variables, i.e. Sample 1 - very small sample but all variables are available for all obs Sample 2 - somewhere in between Sample 3 - all waves and observations, but for many individuals observations on the interesting variables not available. So now a t-test was used to check whether the mean of the dependent variable and some other unaffected key variables differ among these samples. Unfortunately, they do in most cases. I did the t-test in the first instance to prove that for example the dependent variable does not significantly differ between samples. I am not sure what to do with the information, that they do differ. Is there anything I can do?


title: Entropy (individuals equally high on all items vs. high on only few items)

text: I remember something from a stats course many years ago which might be helpful now. I want to distinguish between patients who show symptoms all over the board, vs. patients who have similarly high sum-scores, but due to a few very high symptoms. Symptoms have a scale from 0 to 3. I remember the following procedure (imagine we have 3 symptoms per patient):   1. Calculate the percentage each symptom in relation to the sum-score of that person's symptoms, so s1/sum, s2/sum, s3/sum. Let's call these 3 values p-values.    2. Build natural log of these 3 values. Let's call these p(ln)-values. (If original value is 0, result of this should be 0, or change scale that there are no 0s anymore)   3. Calculate: -2*p*p(ln) per symptom.    4. Sum these up over all symptoms.  Using this procedure, I get the same values if all symptoms are equally high (no matter if they are all 1 or all 3), which is what I want. However, the variance in the value between participants with very equal answering patterns, and only few very high symptoms, are quite minimal, which could be due to the small symptom scale 0-3 (the differences get larger using higher values). I am not sure whether I misremember the name "entropy", or the formula, and would appreciate help. Could I "inflate" the differences, e.g. using 1 10 100 1000 instead of 0 1 2 3 as symptom values?


title: How to interpret coefficients of a multivariate mixed model in lme4 without overall intercept?

text: I'm trying to fit a multivariate (i.e., multiple response) mixed model in `R`. Aside from the `ASReml-r` and `SabreR` packages (which require external software), it seems this is only possible in `MCMCglmm`. In the paper that accompanies the `MCMCglmm` package (pp.6), Jarrod Hadfield describes the process of fitting such a model as like reshaping multiple response variables into one long-format variable and then suppressing the overall intercept. My understanding is that suppressing the intercept changes the interpretation of the coefficient for each level of the response variable to be the mean for that level. Given the above, is it therefore possible to fit a multivariate mixed model using `lme4`? For example:               data(mtcars)     library(reshape2)     mtcars <- melt(mtcars, measure.vars = c("drat", "mpg", "hp"))     library(lme4)     m1 <- lmer(value ~ -1 + variable:gear + variable:carb + (1 | factor(carb)),         data = mtcars)     summary(m1)     #  Linear mixed model fit by REML      #  Formula: value ~ -1 + variable:gear + variable:carb + (1 | factor(carb))      #     Data: mtcars      #   AIC   BIC logLik deviance REMLdev     #   913 933.5 -448.5    920.2     897     #  Random effects:     #   Groups       Name        Variance Std.Dev.     #   factor(carb) (Intercept) 509.89   22.581       #   Residual                 796.21   28.217       #  Number of obs: 96, groups: factor(carb), 6     #       #  Fixed effects:     #                    Estimate Std. Error t value     #  variabledrat:gear  -7.6411     4.4054  -1.734     #  variablempg:gear   -1.2401     4.4054  -0.281     #  variablehp:gear     0.7485     4.4054   0.170     #  variabledrat:carb   5.9783     4.7333   1.263     #  variablempg:carb    3.3779     4.7333   0.714     #  variablehp:carb    43.6594     4.7333   9.224      How would one interpret the coefficients in this model? Would this method also work for generalized linear mixed models?


title: Sign flipping when adding one more variable in regression and with much greater magnitude

text: Basic setup: regression model: $y = \text{constant} +\beta_1x_1+\beta_2x_2+\beta_3x_3+\beta_4x_4+\alpha C+\epsilon$ where C is the vector of control variables. I'm interested in $\beta$ and expect $\beta_1$ and $\beta_2$ to be negative. However, there is multicollinearity problem in the model, the coefficient of correlation is given by, corr($x_1$,$x_2)=$ 0.9345, corr($x_1$,$x_3)=$ 0.1765, corr($x_2$,$x_3)=$ 0.3019. So $x_1$ and $x_2$ are highly correlated, and they should virtually provide the same information. I run three regressions:   1. exclude $x_1$ variable; 2. exclude $x_2$ variable; 3. original model with both $x_1$ and $x_2$.  Results:   For regression 1 and 2, it provide the expected sign for $\beta_2$ and $\beta_1$ respectively and with similar magnitude. And $\beta_2$ and $\beta_1$ are significant in 10% level in both model after I do the HAC correction in standard error. $\beta_3$ is positive but not significant in both model. But for 3, $\beta_1$ has the expected sign, but the sign for $\beta_2$ is positive with the magnitude twice greater than $\beta_1$ in absolute value. And both $\beta_1$ and $\beta_2$ are insignificant. Moreover, the magnitude for $\beta_3$ reduces almost in half compared to regression 1 and 2. My question is: Why in 3, the sign of $\beta_2$ becomes positive and much greater than $\beta_1$ in absolute value? Is there any statistical reason that $\beta_2$ can flip sign and has large magnitude? Or is it because model 1 and 2 suffer omitted variable problem which inflated $\beta_3$ provided $x_2$ has positive effect on y? But then in regression model 1 and 2, both $\beta_2$ and $\beta_1$ should be positive instead of negative, since the total effect of $x_1$ and $x_2$ in regression model 3 is positive.


title: Outcome variable in experimental design

text: I have been asked to try to operationalize our research question methodologically in an experimental design but since this is my first experience doing so I feel unsure and am looking for some help. We are looking at reputation systems (such as that of eBay) where provision of rewards and/or sanctions can stabilize markets and facilitate trade even in the absence of a functioning legal framework. We are specifically interested in analyzing the causal effect of two variables, X and Y on the type of information subjects provide and hypothesize, for example, an interaction effect between X and Y to be related to lower likelihood of sanctioning. I am therefore thinking about using a 2x2 design and my question concerns how to define our outcome variable. A consideration is to let subjects rate their trading partners on a scale from -2 to +2 where they can thus either provide sanctions in terms of negative scores or rewards in terms of positive scores (i.e. -2, -1, 0, +1, +2). Would it be possible to define it in this way even though we are interested in analyzing the effects of X and Y on sanctions (negative ratings -2, -1) separately from their effects on rewarding (positive ratings -1, -2)? And if so, how could we deal with ratings of 0? I would be very grateful to any help or suggestions. Thanks, KML


title: Confidence Interval of estimator for the exponential distribution

text: > We have an exponential distribution > $$f(x)=\frac{e^\frac{-x}{\theta}}{\theta}$$ > > We are told that $n=3$ and that the data are given as $x_1=1, x_2=2.5, > x_3=5.5$ > > a.) Determine an approximate 95% confidence interval based on asymptotic > normality of $\overline{x}$ > > b.) Determine an approximate 95% confidence interval based on the asymptotic > distribution of the deviance $D(\theta)$. The solution does not have to be > exact. Pay attention to the function g, defined as $$g(\theta)= > \frac{3}{\theta}+\ln\theta -\ln3-1$$ My questions: a.) I have to take the sample mean, and then I add up the $z_{ \frac{\alpha}{2} }\frac{\sigma}{\sqrt n}$ thing? With $z_{0.025}=1.96?$ What value do I use for $\sigma$? b.) What deviance $D(\theta)$ mean? The question added a graph, on the x-axis I see theta from 0 to 14 and on the y-axis I see (a probability?) from 0.0 up to 1.0 Has this something to do with the power function of the test? I should have knowledge about Neyman-Pearson, likelihood ratio tests, power function of a test.


title: Simple question Maximum Likelihood

text: I have the following distribution, defined for $0 < x < \theta$, its value is $0$ otherwise. $$f_\theta(x)= \frac{2x}{\theta^2} $$ Find the MLE of $\theta$ I tried: $$\prod_{i=1}^n \frac{2}{\theta^2}x_i =\left(\frac{2}{\theta^2}\right)^n \prod_{i=1}^n x_i $$ Taking the natural logarithm gives us: $$ n \ln{\frac{2}{\theta^2}} + \sum\ln({x_i})=2n \ln\left({\frac{\sqrt2}{\theta}}\right) + \sum\ln({x_i})$$ Taking the derivative with respect to $\theta$: $$ \frac{2n\theta}{\sqrt{2}}=0$$ After this question I get other questions about this MLE (unbiasedness, consistency, sufficiency etc), so I have the feeling that this estimator has to be a 'concrete' value... What's going on? :-)


title: Using glmer, why is my random effect zero?

text: We’ve run a mixed effects logistic regression using the following syntax;               # fit model     fm0 <- glmer(GoalEncoding ~ 1 + Group + (1|Subject) + (1|Item), exp0,                  family = binomial(link="logit"))     # model output     summary(fm0)      Subject and Item are the random effects. We’re getting an odd result which is the coefficient and standard deviation for the subject term are both zero;               Generalized linear mixed model fit by maximum likelihood (Laplace     Approximation) [glmerMod]     Family: binomial  ( logit )     Formula: GoalEncoding ~ 1 + Group + (1 | Subject) + (1 | Item)     Data: exp0          AIC      BIC      logLik deviance df.resid      449.8    465.3   -220.9    441.8      356           Scaled residuals:      Min     1Q Median     3Q    Max      -2.115 -0.785 -0.376  0.805  2.663           Random effects:     Groups  Name        Variance Std.Dev.     Subject (Intercept) 0.000    0.000        Item    (Intercept) 0.801    0.895        Number of obs: 360, groups:  Subject, 30; Item, 12          Fixed effects:                     Estimate Std. Error z value Pr(>|z|)          (Intercept)     -0.0275     0.2843    -0.1     0.92          GroupGeMo.EnMo   1.2060     0.2411     5.0  5.7e-07 ***      ---      Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Correlation of Fixed Effects:                  (Intr)      GroupGM.EnM -0.002      This should not be happening because obviously there is variation across subjects. When we run the same analysis in stata               xtmelogit goal group_num || _all:R.subject || _all:R.item          Note: factor variables specified; option laplace assumed          Refining starting values:           Iteration 0:   log likelihood = -260.60631       Iteration 1:   log likelihood = -252.13724       Iteration 2:   log likelihood = -249.87663            Performing gradient-based optimization:           Iteration 0:   log likelihood = -249.87663       Iteration 1:   log likelihood = -246.38421       Iteration 2:   log likelihood =  -245.2231       Iteration 3:   log likelihood = -240.28537       Iteration 4:   log likelihood = -238.67047       Iteration 5:   log likelihood = -238.65943       Iteration 6:   log likelihood = -238.65942            Mixed-effects logistic regression               Number of obs      =       450     Group variable: _all                            Number of groups   =         1                                                          Obs per group: min =       450                                                                    avg =     450.0                                                                    max =       450          Integration points =   1                        Wald chi2(1)       =     22.62     Log likelihood = -238.65942                     Prob > chi2        =    0.0000          ------------------------------------------------------------------------------             goal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------        group_num |   1.186594    .249484     4.76   0.000     .6976147    1.675574            _cons |  -3.419815   .8008212    -4.27   0.000    -4.989396   -1.850234     ------------------------------------------------------------------------------          ------------------------------------------------------------------------------       Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]     -----------------------------+------------------------------------------------     _all: Identity               |                    sd(R.subject) |   7.18e-07   .3783434             0           .     -----------------------------+------------------------------------------------     _all: Identity               |                      sd(R.trial) |   2.462568   .6226966      1.500201    4.042286     ------------------------------------------------------------------------------     LR test vs. logistic regression:     chi2(2) =   126.75   Prob > chi2 = 0.0000          Note: LR test is conservative and provided only for reference.     Note: log-likelihood calculations are based on the Laplacian approximation.      the results are as expected with a non-zero coefficient / s.e. for the Subject term. Originally we thought this might be something to do with the coding of the Subject term, but changing this from a string to an integer did not make any difference. Obviously the analysis is not working properly, but we are unable to pin down the source of the difficulties. (NB someone else on this forum has been experiencing a similar issue, but this thread remains unanswered link to question) Any comments greatly appreciated!


title: Too good results from EM for gaussian mixture

text: My task is to identify parameters (mean, standard deviation, height) of gaussian peaks in given histogram data with as lowest CV as possible. Number of peaks and approximate means are known (pointed by user). Standard approach uses regression to fit gaussian mixture to the data. Results are listed below: ![enter image description here](http://i.stack.imgur.com/Ejy4H.png) ![enter image description here](http://i.stack.imgur.com/QQJBj.png) I've tryed to employ EM algorithm in the following way:   * At first, I didn't care about heights of peaks (used normalized values) and I was searching only for means and standard deviations using EM.   * Second, I fixed means and standard deviations found by EM and used regression to find height of each gaussian.   * I received following results: ![enter image description here](http://i.stack.imgur.com/oqK6u.png) ![enter image description here](http://i.stack.imgur.com/1DAEi.png) I want to ask:   * If my approach is theoretically correct.   * If those results from second approach could be (by sense) valid. Except the last peak, of cause, which had so few points.  Chi squared result is logically bad for second solution; however, the data could contain outliers. I would be very happy if those CV values were correct, but I am not sure, if I didn't do something wrong. **UPDATE:** I try to explain briefly, why I am so surprised... These measured data were collected from cytometer and correspond to particle (biological cell) count (y axis) per some fluorescence intensity (x axis).. This intensity correlates with biological cell size, and the goal is to measure cell size of each examined type of plant (one peak correspond to one type of plant). It is expected that not all cells of same plant share same size, but we know that these sizes (of one plant) are normally distributed. During measurement process many outliers and false signals are collected and mixed with the original data (we know that the real CV should be much lower), so we want to reach as lowest CV as possible to be sure about cell-type size. Cytometric softwares that we tested uses the method that i described first (non-linear least squares regression). I've tryed EM, received low CVs, however received a huge difference between my model and data (can be seen in image or ChiSqr). Now I am wondering, if I didn't missed something, because if my approach was correct, then common cytometric softwares would probably use it.


title: Why does the log-likelihood change when a variable is linearly transformed in a hierarchical model?

text: I ran into (what I think is) an inconsistency when running a random-intercept model (using the _lmer_ function in the _lme4_ package in R). Here is what I do: I first run a model with a set of covariates; then I run the same model re-scaling (linearly transforming) one of the regressors. To my knowledge, this should change **only the coefficient** of the variable that is linearly transformed. And indeed, this is what happens when I run this "experiment" with a simple linear regression model and with a logistic model. This code replicates the "normal" behaviour:               # Create three random independent variables     x1 <- rnorm(20)     x2 <- rnorm(20)     x3 <- as.factor(sample(0:2, 20, replace = TRUE))     # Their random coefficients     coef1 <- runif(1, -1, 1)     coef2 <- runif(1, -1, 1)     # Create a continuous dependent variable and a binomial one     y1 <- coef1 * x1 + coef2 * x2 + runif(20)     y2 <- y1     y2[which(y1 > quantile(y1, 0.5))] <- 1     y2[which(y1 <= quantile(y1, 0.5))] <- 0     # Finally, a linear transformation of x1     x1.trans <- x1*3      So, let us run an OLS model:               lm <- lm(y1 ~ x1 + x2 + x3)     summary(lm)     # OLS model with one variable linearly transformed     lm.bis <- lm(y1 ~ x1.trans + x2 + x3)     summary(lm.bis)      The coefficients of _x1_ and _x1.trans_ are different, **but the R-square of the two models is the same** :               summary(lm)$r.sq == summary(lm.bis)$r.sq      The same with a logistic model:               logm <- glm(y2 ~ x1 + x2, family="binomial")     summary(logm)     logm.bis <- glm(y2 ~ x1.trans + x2, family="binomial")     summary(logm.bis)      Even in this case, **the log-likelihood of the two models is the same** :               logLik(logm) == logLik(logm.bis)      So far, so good. However, when I do the same with a hierarchical model, **the log-likelihood (and consequently the AIC and BIC) of the two models are different** , although the coefficient of the transformed variable remains significant with the same z value and the other coefficients are the same.               # Multilevel model     mm <- lmer(y1 ~ x1 + x2 + (1 | x3))     summary(mm)     mm.bis <- lmer(y1 ~ x1.trans + x2 + (1 | x3))     summary(mm.bis)     logLik(mm) == logLik(mm.bis) ### FALSE! ###      Why? Also the "REML criterion at convergence" is obviously different. I don't understand this result. This is probably due to my moderate knowledge of the math of hierarchical models. I'd be very happy if some of you could show me what's the trick here. Since we then use AIC and BIC to compare models, I am puzzled by the fact that a simple transformation that shouldn't change anything makes one model better (or worse) than another.


title: Yet another "Bayesian vs Maximum Likelihood" question

text: In the fully Bayesian approach, the predictive distribution is: $$ P( Y|X ) = \int P(\theta | X ) P( Y | \theta ) d\theta $$ When the integral is difficult to compute, we might resort to the Maximum Likelihood approach, and approximate the predictive distribution as follows: $$ P( Y|X ) \approx P( \hat \theta_{ML} | X ) P( Y | \hat \theta_{ML} ) $$ where $ \theta_{ML} $ is the MLE of the parameter $\theta$. When you look at this, the ML approach looks like an awfully bad approximation. Why does it work not so badly in practice?


title: Derivation of Zhou (1996) volatility estimator

text: Does anyone know how to derive the Variance of Bin Zhou's volatility estimator (Theorem 1) in 'High-Frequency Data and Volatility in Foreign-Exchange Rates' (1996)? PDF Link attached. Any help would be greatly appreciated.


title: how to initialise this distribution

text: I am reading a paper by Tom Minka and trying to follow the steps to perform the updates for the Expectation Propagation algorithm. I am having trouble getting the initialisation of the posterior distribution going. So, the first step involves writing out your joint distribution as a product of factors. In my case, I can do the following: $$ P(\theta, D) = g(\theta) \prod_i f_i(\theta) $$ where $g(\theta)$ is a multivariate normal with a certain mean ($\mu)$ and a covariance matrix ($\Sigma$). Each of the factors $f_i$ can be written as a normal distribution with a mean ($m$) and a standard deviation ($\sigma$). Now, the paper (here) says as the first step, compute the approximate posterior as: $$ q(\theta) = \frac{g(\theta) \prod_i f_i(\theta)}{\int g(\theta) \prod_i f_i(\theta) d\theta} $$ Now, I am really confused as how I can do this operation! So, I can initialise the $g(\theta)$ and each of the $f_i(\theta)$ terms with some initial values but how do I compute $q(\theta)$ in the right form. More specifically, if I assume $q(\theta)$ will take the form of a multivariate Gaussian. How do I compute its mean and covariance from this expression? I have been looking online and it seems that there are some Gaussian identities one can use but none of the papers actually make it explicit as to which ones to use.


title: R - Network Connectivity over the time

text: I have to analyze if a network (represented as a graph) mantains its connectivity over the time. The data that I've obtained from the simulation is given here, where L1 points out the vertex of the graph, and the column `value` stands for the degree of each node. This data is moreover obtained from this step of the computation in R:               m <- melt(x,id="time")      The plotting of               ggplot(m,aes(x=time,y=value)) +stat_summary(fun.data = "mean_sdl", geom = "smooth")+stat_summary(fun.data = "mean_sdl", geom = "errorbar")      is given in this following picture: ![enter image description here](http://i.stack.imgur.com/M8axH.png) Seeing that the result had many errors and that I wasn't able to check where the connectivity remained stable or not, I've tried to evaluate for each node's degree its `acf` over the time (as a temporal series) and hence I've obtained the other following result with the same plotting command: ![enter image description here](http://i.stack.imgur.com/HnaZ1.png) This didn't help to lower the error. I'm asking if there's a tecnique to get lower errors and by which determine if a network keeps the same connectivity over the time. Could a frequency analysis over the number of nodes with a same degree over the time helpful and lower my errors? Thanks in advance.


title: Range of a centered, reduced sample drawn from a bounded domain

text: Assume you draw $N$ values from a bounded set, say $[0, 1]$ without loss of generality. If I am right, the centered values $X-m$ will remain in range $[{1\over N}-1, 1-{1\over N}]$ (by considering the most imbalanced 0/1 cases). And the standard deviation in range $[0,{1 \over 2}]$ (even $N$) or $[0,{1\over 2}\sqrt{1-{1\over N^2}}]$ (odd $N$) (by considering the most balanced 0/1 cases). But what is the range of the centered-reduced values ${X-m\over s}$ ? How can I address this ?


title: Making sense of the allele frequency weighted genetic distance

text: I found one formula to calculate the pairwise distance between samples according to their SNPs in the following paper: Siu, Jin and Xiong, Manifold Learning for Human Population Structure Studies, PLoS One (2012) 7(1), p. 13. The formula is $d = 1-\frac{1}{2k}\sum_{j=1}^k (\frac{1}{P_{A_j}}IBS_{A_j}+\frac{1}{P_{a_j}}IBS_{a_j}+\frac{1}{P_{N_j}}IBS_{N_j})$ where $IBS_{A_j},IBS_{a_j},IBS_{N_j}$ be the number of copies of alleles $A_j,a_j$ and missing allele $N_j$ at the $j$-th SNP shared by the pair of individuals, respectively, and $P_{A_j},P_{a_j},P_{N_j}$ are the frequencies of the alleles $A_j,a_j,$ and $N_j$, respectively, and $k$ is the number of SNPs. For simplicity, let's assume that there is no missing values and $k=1$. We would have cases like AA-AA, Aa-Aa, aa-aa, AA-Aa, Aa-aa, AA-aa. The distance for the first three case should be $0$. Let's use the formula to calculate the case for Aa-Aa: we should have $IBS_A=1,IBS_a=1,P_A=\frac{1}{2},P_a=\frac{1}{2}$, then the formula gives me a negative distance. How could that be possible? To my understanding, $d$ can equal to $0,\frac{2}{3},1$. We have $\frac{2}{3}$ because 'weighted'? Is there any other formula to calculate the pairwise distance according to sample's SNPs?


title: Multiple Variable Ranking

text: Before i start, i should mention i am pursuing my under-grad degree, and hence i have little/basic knowledge about machine learning and statistics. I've been meaning to make an application which takes a 5-point rating of a particular item based on a few variables. > **Example:** > > **Some Food Item** (Gets an average rating based on the rating of its > variables - out of 5) > > Variables:(Each variable is rated on a 5-point score): > > Presentation, Taste, Value for Money, etc. This is all the User Input ill have. I am searching for a Rating System/Algorithm that will help me take into consideration every variable under each item. **Research**   * Bayesian System. **The Problem** : It doesnt take into to consideration the multiple variables for rating.   * Good Grief Algorithm. **The Problem** : I didn't understand how it is implemented, in my scenario here. **The Question:**   1. Is Bayesian System the right way? If yes, how do I consider the multiple variables?   2. The (external link) PDF of Good Grief Algorithm, has the exact problem statement as mine, but i do not understand the implementation. If it's the right way, can you break down the implementation, in simple pseudo code(or math formulas)?   3. Is there is a better/simpler method of rating, please explain it in the answers. P.S. Again, i'd say, I'm not exactly a math genius, so it would be better if the answers are elaborate, and broken down to the most basic formula. Pseudo Code is preferred. Thank you.


title: do i need to add a variable to my regression?

text: I am doing multiple regression on a subset of the NLYS79 dataset, namely a subset containing 540 respondents. And am interested in the significance of race upon earnings. My variables here are: years in school, ASVABC(a measure of intelligence), ethnicity is black, ethnicity is white, female, tenure, hours per week worked and married*male. $$\log(\text{earnings}) = \alpha +\beta_1\text{ASVABC}+\beta_2\text{ethblack}+ \beta_3\text{ethwhite}+\beta_4\text{female} \\\\+\beta_5\text{Schooling}+\beta_6 \text{tenure}+\beta_7\text{hours}+\beta_8\text{married*male}$$ For all these variables we have a source containing information about how this affects the wage. Now my question is: all my variables except ethblack/ethwhite are significant (|t-statistic| > 2) am i safe to conclude anything from my data? How do I know that there is something wrong with my regression equation? In essence my result is something I want: which I have now, but the whole "ommission of variable" thing is puzzling me at the moment. I know there are things wrong with my dataset, for example I do not have a 'fair' distribution of the races black, white and hispanic are 63, 599, 34 respectively.


title: What is the difference between using the multiplication rule or using Venn diagram subtraction for probability?

text: Take two problems:   1. Andrew is 35, and the probability he will be alive in 10 years is .72. Ellen is 35, for her, .92. Assuming these are independent, what is the probability they both will be alive in 10 years?  Answer: .66, Method: Use the multiplication rule .72*.92 = .66   2. Suppose your street has two traffic lights. The chance that the first light is red is .40, and the second light .30. The chance of them being red at the same time is .10. What is the the probability that neither light is red? Answer: .40, Method: Subtract .1 from .4 and .3, and then subtract all 3 from 1.0 = .4 **Why doesn't the multiplication rule work for #2?** In other words, why can't I multiply .40 and .30 to get .12? And then furthermore, if .12 is the chance they will both be red at the same time shouldn't (1-.12)=.88 be the chance that neither is on?


title: Obtaining a log-normal waiting time via sequential exponential or gamma distributions - is it possible?

text: Assume you've use a accelerated failure time model to find that the transition of subjects from State A to State B is log-normally distributed with parameters $\mu$ = X and $\sigma$ = Y. This now needs to be used in a differential equation model as the rate at which subjects move from State A to State B. However, just using $\mu$ as the average probability of moving from A to B results in an exponentially distributed waiting time, not a log-normally distributed time. I know you can however use a series of sequential independent exponential distributions to obtain an overall gamma distributed waiting time. For example, if it takes 2 days to move from A to B, four equally spaced exponential distributions results in an overall gamma distributed waiting time with $\kappa$ = 4 and $\theta$ = 2. The question is, is there a gamma distribution whose shape and scale parameters approximate a log-normal? I know if I used many exponential distributions for a high $\kappa$ the central limit theorem allows the gamma distribution to approximate a log-normal, but I'm not sure if there's a way to obtain a log-normal waiting time from that. Essentially, is there some $\kappa$ and some $\theta$ that yields something approximating a log-normal with $\mu$ = X and $\sigma$ = Y?


title: Pooling the results of random hot-deck imputation

text: I am using random hot-deck imputation on a repeated measures dataset. I am tempted to use Rubin's rules for pooling the results of multiple imputation, in particular for regression coefficients. Intuitively it seems the average of the coefficient estimates could be used, but I really don't have any insight about pooling the standard errors and I have not seen any literature about this (Little & Rubin's book seems to be silent on the matter, unless I have missed something) How can I pool regression coefficient estimates and their standard errors when performing random hot-deck imputation ? Pragmatic advice, theoretical justification and/or references to the literature would be most welcome. **Edit** : To clarify, by "random hot deck imputation", I mean: > Hot deck imputation involves replacing missing values of one or more > variables for a non-respondent (called the recipient) with observed values > from a respondent (the donor) that is similar to the non-respondent with > respect to characteristics observed by both cases. In some versions, the > donor is selected randomly from a set of potential donors, which we call the > donor pool; we call these methods random hot deck methods http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3130338/


title: Imputation for a zero-inflated negative binomial mixed effects model

text: I am working with a dataset of repeated (x4) observations on 100 subjects. The outcome is zero-inflated and the data appears to be modelled well by a mixed effects zero-inflated negative binomial model with random intercepts for subjects. However, approximately 20% of the outcome variable is missing, so I am investigating imputation methods that can deal with this. There are 2 covariates in the analysis model and a further 6 auxiliary variables that can be used for imputation - all of which have very low levels of missingness. I have tried the `mice` package in R which supports random effects, but does not support the negative binomial distribution. The only method I have found useful so far, is random hot-deck imputation, which I have tried in R (with the package `StatMatch`) and Stata (with the package `hotdeck`) both of which seem to produce reasonable results Is there any other package/system that can be used to impute values for this kind of model ? I have experience with R and Stata, and SAS is also available to me (though I don't have any experience with it)


title: Variation explained in ordinal logistic regression models

text: I have made these three ordinal logistic regression models:               model1 <-  polr(as.factor(carb) ~ mpg,  Hess = T, mtcars)     model2 <-  polr(as.factor(carb) ~ hp,   Hess = T, mtcars)     model3 <-  polr(as.factor(carb) ~ drat, Hess = T, mtcars)      To figure out if the models are a good fit to the data, I calculated the proportion of variation explained like this:               model_null <-  polr(as.factor(carb) ~ 1, Hess = T, mtcars)          1-(model1$deviance/model_null$deviance)     0.1512784     1-(model2$deviance/model_null$deviance)     0.2520109     1-(model3$deviance/model_null$deviance)     0.003453936      Questions:   1. Why doesn't summary give null deviance?    2. Have I calculated the proportion of variation explained correctly?   3. Am I right in saying `model1` and `model3` explain little variation in `carb`, but `model2` explains 25% of the variation in `carb`?


title: Comparing two Bayesian models

text: To compare two Bayesian models for the same data, I plotted the fitted values of the first model versus the fitted values of the second model. The plot indicates a perfect similarity. Relying on this plot, I said there is no significant difference between these two models. However, the difference between the Deviance information criterion (DIC) values of each model is 85, which I think shows that there is a significant difference between these two models. So, the plot and the DIC, I suppose, are in contrast. What does everyone else think?


title: string clustering: similarity criterion

text: I have a set of strings of dimension $10,000$. I want to group similar strings together in one group, perform clustering. As string metric, I am using the `Levenshtein distance`. * * * Simply, with the `Levenshtein distance` I'll just compute distance between $2$ strings and then by using a threshold the clustering algorithm will make the decision if they can be grouped or not. This is not enough. I am looking for a special measure to study the relation between the strings. * * * **For example:** **door** and **entrance** wont be grouped together if I just compute the `Levenshtein distance`, in fact there is nothing in common between these $2$ words. But logically they are connected and can be grouped together, since **door** and **entrance** are basically the same. * * *   * Have you ever come across such a problem?   * Is it similar to the **Semantic similarity** measure?


title: Clustered standard errors and multi-level models

text: Stata allows estimating clustered standard errors in models with fixed effects but not in models random effects? Why is this? By clustered standard errors, I mean clustering as done by stata's cluster command (and as advocated in Bertrand, Duflo and Mullainathan). By fixed effects and random effects, I mean varying-intercept. I have not considered varying slope.


title: How should I measure the relationship between a variable containing text and a binomial variable?

text: I have data in this format:               Words,Source,Result     pestology Gomel, cheap, 0     cocreating, cheap, 1     munitioner impersonating nonextinct, cheap, 1     Kolomna, expensive, 1     Enyo's snakemouth, expensive, 0     blueberries backare farriers, cheap, 0     markets rafales, cheap, 0     ...      http://bit.ly/h2ynoG How should I determine whether the phrases in column one can explain the Result value? What if I want to determine which individual words in the phrase are correlated with successful results? Also, how should I determine whether the data in which Source is cheap is representative of all of the data? Thanks.


title: No lrtest() function available in lmtest

text: I have got a strange problem with using functions which are supposed to be in the libraries lmtest & language R (downloaded from cran.r-project recently), namely lrtest() and pvals.fnc() - however, they are not available. I have upgraded my R-version to 2.12.2 but it has not solved the problem.               library(lmtest)     lrtest()     returns: Error: could not find function "lrtest"      Any help would be greatly appreciated. T


title: How can I test the quality of a RNG?

text: On electronics.stackexchange we had a question about constructing a True Random Number Generator. Since the method relies on noise, which isn't deterministic, the only way to test the quality of the RNG seems to be empirical.   Not being a statistician I suggested to test a long bit sequence for normality, but I have no idea what outcome is acceptable and what isn't. For instance counting single bits I guess we have a thumbs-up for a 499500/500500 distribution, but when is the ratio too much skewed to be acceptable? Same for 2-bit and longer sequences.   Of course, if testing for normality is a Bad Idea™, I'd like to hear it, including better alternatives. **edit**   Diehard was mentioned a few times, but I'm not sure this answers my question. The normality test should give a normal distribution, with better bell curve approximations for longer sequences. Diehard seems also to have tests which either should result in normal or exponential distributions. But my question remains: how do I judge the results? Just by looking at the curve and discern a bell curve in it? To get back at my first example, a 499500/500500 distribution is definitely OK, and a 950000/50000 distribution definitely is a no-no, so where does the switch-over happen?


title: Looking for stats/probability practice problems with data and solutions

text: I am looking for a self-study site on the web that will allow me to verify my understanding of some basic probability and stats concepts and operations. What I would like is a site with data and problems (along with the solutions) that I could use to practice my skills (such as regression methods/analysis, pairwise t-tests, computation and interpretation of confidence intervals, etc.). I've seen workbooks like this, but the data is not available for me to import into R or Excel. Surely there must be a site like this on the internet, arranged by topics for study, but somehow I have not been able to find it. Also, if anyone knows of a good workbook that comes with data (on CD or via web) I'd be interested in it too. Re tools, I'll be using R and Excel - so pointers re these would be great too. Thanks!


title: Assessing the effect of adding a variable using stepwise forward logistic regression using Stata?

text: I'd really appreciate help using Stata to perform a manual stepwise forward logistic regression. I have 37 biologically plausible, statistically significant categorical variables linked to disease outcome. I need to end up with a final multivariable model. I've added the first variable (most significant/most plausible) with corresponding OR output.               xi:logistic outcome i.variable1          . xi:logistic casecontrol i.breed_groupall     i.breed_group~l   _Ibreed_gro_0-7     (naturally coded; _Ibreed_gro_0 omitted)          Logistic regression                               Number of obs   =        995                                                       LR chi2(6)      =      83.87                                                       Prob > chi2     =     0.0000     Log likelihood = -422.36813                       Pseudo R2       =     0.0903          ------------------------------------------------------------------------------      casecontrol | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------     _Ibreed_gr~2 |   1.757143   .6861797     1.44   0.149      .817345    3.777537     _Ibreed_gr~3 |   1.952381   .8811439     1.48   0.138     .8061249    4.728537     _Ibreed_gr~4 |   1.464286    .530121     1.05   0.292     .7202148    2.977074     _Ibreed_gr~5 |   6.192708   1.779109     6.35   0.000     3.526453    10.87485     _Ibreed_gr~6 |   3.880357   1.103611     4.77   0.000     2.222193    6.775816     _Ibreed_gr~7 |    .636646   .2555236    -1.13   0.261     .2899083    1.398091     ------------------------------------------------------------------------------      What do I look for to see if adding the second variable I choose means that both variables should stay in, when, for example I type;               xi:logistic outcome i.variable1 i.variable2               . xi:logistic casecontrol i.breed_groupall i.height_category     i.breed_group~l   _Ibreed_gro_0-7     (naturally coded; _Ibreed_gro_0 omitted)     i.height_cate~y   _Iheight_ca_0-4     (naturally coded; _Iheight_ca_0 omitted)          Logistic regression                               Number of obs   =        992                                                       LR chi2(10)     =     132.25                                                       Prob > chi2     =     0.0000     Log likelihood = -396.05629                       Pseudo R2       =     0.1431          ------------------------------------------------------------------------------      casecontrol | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------     _Ibreed_gr~2 |   1.002262   .4145417     0.01   0.996     .4455736    2.254465     _Ibreed_gr~3 |   1.185087    .557553     0.36   0.718     .4712828    2.980017     _Ibreed_gr~4 |   1.774162   .6616502     1.54   0.124     .8541793    3.685001     _Ibreed_gr~5 |   1.452416    .541494     1.00   0.317     .6994292     3.01605     _Ibreed_gr~6 |   1.256098   .4377793     0.65   0.513     .6343958    2.487064     _Ibreed_gr~7 |   .4238999   .1777699    -2.05   0.041     .1863361    .9643388     _Iheight_c~1 |   2.780857    .918967     3.09   0.002     1.455088    5.314572     _Iheight_c~2 |   5.402833   2.246817     4.06   0.000     2.391342    12.20679     _Iheight_c~3 |   13.50715   5.989787     5.87   0.000     5.663642    32.21303     _Iheight_c~4 |   16.85605   8.674745     5.49   0.000     6.147467    46.21846     ------------------------------------------------------------------------------      How do I know if I want to keep one, or both of these variables, or that one, or both of them is no use to me? If for example, I want to keep both of these and add the 3rd variable, how do I know which?               . xi:logistic casecontrol i.height_category i.breed_groupall i.combinedweight          i.height_cate~y   _Iheight_ca_0-4     (naturally coded; _Iheight_ca_0 omitted)     i.breed_group~l   _Ibreed_gro_0-7     (naturally coded; _Ibreed_gro_0 omitted)     i.combinedwei~t   _Icombinedw_0-5     (naturally coded; _Icombinedw_0 omitted)          Logistic regression                               Number of obs   =        891                                                       LR chi2(14)     =     123.58                                                       Prob > chi2     =     0.0000     Log likelihood = -346.82026                       Pseudo R2       =     0.1512          ------------------------------------------------------------------------------      casecontrol | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------     _Iheight_c~1 |   3.397418   1.235922     3.36   0.001     1.665316    6.931087     _Iheight_c~2 |   6.321891   3.119652     3.74   0.000     2.403289    16.62983     _Iheight_c~3 |   14.36312   7.797083     4.91   0.000     4.956448    41.62242     _Iheight_c~4 |   23.39157   15.34571     4.81   0.000      6.46607    84.62101     _Ibreed_gr~2 |   .7339708   .3482777    -0.65   0.515     .2895834    1.860303     _Ibreed_gr~3 |   1.060443   .5318773     0.12   0.907      .396787    2.834113     _Ibreed_gr~4 |   1.644423   .6502556     1.26   0.208      .757569    3.569479     _Ibreed_gr~5 |   1.246412   .4940969     0.56   0.578     .5731024     2.71076     _Ibreed_gr~6 |   1.262449   .4661338     0.63   0.528     .6122442    2.603172     _Ibreed_gr~7 |    .401331   .1782529    -2.06   0.040     .1680497    .9584456     _Icombined~1 |    1.21903   .7160999     0.34   0.736     .3854692    3.855132     _Icombined~2 |   1.238685   .3967314     0.67   0.504     .6612025    2.320531     _Icombined~4 |   1.764532   .5970935     1.68   0.093     .9090641    3.425031     _Icombined~5 |   2.107871   1.118772     1.40   0.160     .7448366    5.965227     ------------------------------------------------------------------------------      Any help would be gratefully received. (I am using STATA v9.1 I believe.)


title: Relationship between MDL and "difficulty of learning from data"

text: While trying to make sense of MDL and stochastic complexity, I found this previous question: Measures of model complexity, in which Yaroslav Bulatov defines model complexity as "how hard it is to learn from limited data." It is not clear to me how Minimal Description Length (MDL) measures this. What I am looking for is some sort of probability inequality (analagous to the VC upper bound) which relates the "code length" of a model with its worst case behavior on fitting data generated by itself. If such a concrete result cannot be found in the literature, even an empirical example would be enlightening.


title: Measuring dispersion of tokens within a text file

text: I am working on a code analyzer application. It is essentially a piece of software that parses and interprets other programs' code and comes up with various metrics, findings, statistics, and ultimately draws performance warnings and re-engineering proposals. It is still in its infancy and what I am currently doing is parsing source code files for variable declaration, assignment, and references. That part works fine, I am able to accurately extract all variable occurrences, their positions and identify them in one of the three aforementioned categories. My next step is to identify and measure clustering of each variable within the code with the intent of identifying distinct functional areas that can be segregated into its own distinct programmatic block (i.e. a function but not necessarily). So it is kind of a simple attempt at `machine learning`. The data at my disposal consists of the file length and a hierarchical tree of variables, line numbers at which they occur, and how many times they occur on each line. E.g               FOO->5 = 1     FOO->17 = 2     FOO->32 = 1     BAR->6 = 1     BAR->55 = 1     etc.      meaning that `FOO` occurs on lines 5 and 32 once, 17 twice on line 2 (example of a red flag for code smell), `BAR` once on 6 and 55 each etc. What I want to do is measure dispersion of each variable but across the whole file. The fact it needs to be measured against the whole file is what made me think the simple `standard deviation` won't do because that measures only the dispersion of distributions within itself, i.e. there is another parameter involved. So if the source code file length is say 500 lines, and there are half a dozen variables scattered in the code block 125-200, with each variable appearing seldom or not at all outside that range, I want my metric to capture that and identify the code block 125-200 as a good candidate for separation into its own unit. Conversely, if a variable is uniformly scattered across the whole file and is kind of ubiquitous, I want to recognize that and designate it as a global variable. My question is: is there a standard statistical formula/measure that would represent the degree of dispersion of the distribution of occurrences across an area (in this case variables in a file), something like a value-added standard deviation?


title: General definition of stochastic processes

text: I'm trying to understand the basic concept of random processes.I already understood that a continuous-time random process is defined by X(¥,t) where ¥ is each element in the Sample Space and t is the time index. By this definition, we can see the random process as an ensemble of sample realizations ( time-functions defined for each outcome in the Sample Space) or as a set of random variables ( each one defined for a specific time t ). In my mind, the value of the time-function ( sample realizations ) X(¥1,t) reflects the changing value of the probability of outcome ¥1 being choosen in the sample space, is that correct ? For example, we have a boy B in the sample space representing outcome ¥1.Having this outcome ¥1 at time t=0 equals 30%.As the time passes, the sample spaces gets more crowded and at time t=5, the probabily of choosing boy B equals 20%. Can i think that the time-function X(¥1,t) will equal 0,3 for t=0 and 0,2 for t=5 ? Does X(¥1,0)=0,3 and X(¥1,5)=0,2 ? Thanks in advance


title: What is the name of (and alternatives to) this Bayesian point-estimate?

text: Assume that we have a Markovian environment that generates at every time step an event $A$ with probability $p^*$ and an event $B$ otherwise. Now suppose you are a Bayesian agent that wants to learn $p^*$ from the observations you've made so far. If we let our hypothesis class be $\\{H_x\\}_{0 \leq x \leq 1}$ and our initial prior be the uniform distribution ($f_0(x) = 1$), then according to Bayes rule we will update at each time step as: $$ f_{t+1} = \begin{cases} \frac{x f_t(x)}{\int_0^1 x f_t(x) dx} & \text{if we saw A} \\\ \frac{(1 - x) f_t(x)}{\int_0^1 (1 - x) f_t(x) dx} & \text{if we saw B} \end{cases} $$ The best estimate (at time $t$) $p_t$ for $p^*$ will then be the mean of $f_t$. I can show that tracking $p_t$ reduces to counting, in the sense that if by time $t$ you saw $a$ many $A$ events and $b = t - a$ many $B$ events then $$p_t = \frac{a + 1}{(a + 1) + (b + 1)}$$ This seems like a basic result I would have seen in an introductory statistics course (which I unfortunately do not have). **What is this called?** If an agent was using the above procedure to estimate $p^*$ what kind of inference, regression, or estimation would they be doing? Further, to make the answer not just a one-liner. **What are some alternatives?** In particular, what kind of assumptions would I need to make to have $p_t = \frac{a}{a + b}$? I.e. how can I eliminate the uniformity bias caused by my choice of prior?


title: hypothesis testing using poisson distribution

text: At a nuclear plant great care is taken to measure the employees health.These are the number of visits made by each of the 10 employees to the doctor during a calender year. 3,6,5,7,4,2,3,5,1,4 Assuming the number of visits made by employee has a poisson distribution ,test the hypothesis that the annual mean per employee is greater than 3. I am using the graphical method and I am not sure of which p[X=x] should i consider. X: no.of visits by each employee to the doctor.               H0:lambda=3     H1:lambda>3     X follows a Poisson(3)      Then what is the probability that I should check? What I did was as the average of sample data is 4.73636. Therfore calculated p[X>=4] and checked if it was in the critical region. Is this the correct probability to calculate? In a poisson distribution the expected value is calculated as x*p[X=x] right?Not as (sigma x*f(x))/(sigma x)


title: Hypothesis testing with exponential family

text: I'm interested in running hypothesis tests for a variety of members of the exponential family with continuous support, for different values of the parameter/s, for a sample of n i.i.d random variables (distributed according to a particular member and parametrization of the exponential family). In general, it seems the sufficient statistic (at least for a single parameter) is a sum of the variables in the sample. I know this can be obtained at least theoretically from n-fold convolution, but in general, there does not seem to be a nice form for the distribution of this sum, which I will need to compute p-values. so I'm wondering what distributions have known distributions for the sum of n i.i.d variables, and if there is a good reference for doing hypothesis testing with exponential family. I know that:   * sums of iid exponentials are distributed Erlang   * sums of iid normals are distributed normal   * sums of iid gamma are distributed gamma (but what about product of gamma, which is the sufficient statistic for the shape parameter?) So what about beta, pareto, log normal, and so on? Is the difficulty of this the reason why no one does hypothesis tests anymore?


title: False convergence warning message in lmer

text: I will try to explain my data as good as possible. So we taged 13 different whales with a tag that records time, depth, speed, angle of descent and ascent 25 samples every second. The normal diving behaviour of these animals is one deep dive of one hour to 1200 meters followed by a series of 3-7 shallow dives of 20 minutes up to 300 m. Because the tag not always stays the same time in each animal, my data is unbalanced and some tag records have one deep dive and 6 shallow dives while other records have 7 deep dives and 26 shallow dives. I divided each dive in units of 30 seconds. for each unit I have the next data: whale number, dive number, total number of fluke strokes in the 30 seconds unit of analysis, mean of the sin of the angle during the 30 seconds unit, swim speed, dive type(ascent or descent), dive direction( if it is a descent or an ascent) and time since the start of the dive and finally my variable response which is presence or absence of one type of fluke stroke called stroke type B. I think there has to be some autocorrelation between each 30 seconds unit of analyis and need to include it in my model but do not know how! I am interested to know what affects the presence or absence of the type B stroke (which is a binomial variable with 0 and 1) so I decide to use a binomial glmm with whale number as a random effect. I included as well dive number within whale as a random effect. here is the model.               glmm114<-lmer (StrokeB~ Time * Depth+SINP+flukes*Depth+speed+(1|whale_number)+(0+dive_number|whale_number),data=Luciadeepas, family = binomial)      this is my final model after taking out the non significative variables, the problem is that due to the interaction a problem appears saying               The false convergence warning message (8)      I looked on internet and it says is a common problem, and some people says that it doesnt make any change in the output while others says that each variable has to be divided by 100. but when I do this then the variables that become significant doesnt make any sense.