title: Interpreting parameter estimates when multicategorical predictors

text: I would like some advice on how to interpret parameter estimates with a continuous dependent variable (healthy food consumption) and a 6-category predictor (income, from 1 = lowest income to   6 = highest income band in GLM. My output is as follows: Parameter Estimates   Dependent Variable: Healthy $$\begin{array}{c|cccccc}&B&Std. Error&t&Sig.&Lower Bound&Upper Bound\\\hline Intercept& 7.578& 772&9.813&.000&6.059&9.097\\ [income=1]&-1.633 & .826&-1.977&.049&-3.259&-.008\\ [income=2]&-1.310 & .801&-1.635&.103&-2.886&.266\\ [income=3]&-.468 & .831& -.563&.574&-2.102&1.166\\ [income=4]&.664& .927 & .716&.474&-1.159&2.486\\ [income=5]&-.906 & .980& -.924&.356&-2.834&1.022\\ [income=6]&0^\rm a\end{array}$$$\rm a$. This parameter is set to zero because it is redundant. If I plot income in a linear regression model as continuous, I get the following, but I'm not sure it's correct to do that: unstandardised coeff. (income) = .415, sig. = .000. Which would be the best way to intrepret the results? For the first model (GLM), the only significant result I get is for income = 1, which is lowest income, and the B slope is -1.633. I'm not sure how to interpret that. If I would look at the second model (but not sure if it's correct to plot the categorical predictor like that in linear regression), it would suggest the higher the income is the more healthy food people eat. Thank you!


title: Redundant Training Samples

text: I am building a prediction model based on SVM model with RBF kernel. The training samples represent features obtained from genome regions (called binding sites) that are targeted (regulated) by proteins. The model will learn to distinguish between regions targeted by proteins and regions that are not. (I created positive and negative data). My question is that the positive data contains redundant training samples (i.e. identical binding sites at different locations of the genome) so, shall I keep the redundant training samples (to make the model learn more about a specific pattern of targeted region) or I should keep only the unique training samples? For example: I have in the training samples 10,000 samples represent 10,000 different target regions (binding sites) of proteins. If I take the unique, it will be around 7,000 samples. (Because the binding sites might be identical in different locations)


title: Seasonal component for a series that is only one and a half years

text: I have a times series for temperature and I am trying to decompose it into a seasonal + trend + stochastic process to look at what model would fit for the stochastic process. However, I only have data for 333 days, and I assume my data should include a periodic phenomena, with the period being one year, or 365 days. How can I eliminate the trend and seasonal component when I don't have data for the entire assumed period? Is there a way to implement the method in R? Here is the series http://www.wikiupload.com/XD5FHJFC7S857MG EDIT: Plot of this data by @Penguin_Knight ![enter image description here](http://i.stack.imgur.com/j964H.png)


title: Estimating/predicting maximum of sum process for modeling number of cars on the road?

text: I have two series of information about (for example) the amount of cars on a two roads over a 1 mile stretch. I have collected 4 outcomes: Number of cars passing through an intersection at 1, 10, and 60 minutes as well as greatest number of cars in the area in the last hour. I'd like to model the maximum of the sum of the number of cars on both roads in an hour for a given hour. Adding an example for clarity (as per comments): If (over a 20 minute period, sampled every minutes), my observations are `2,4,7,4,8` and `5,4,3,2,1`, the maximum over both roads is 10 (7 + 3), not 13 (5 + 8), because the queue specific maximums happened at different times.


title: Prove that the variance of $b$ is as follows

text: This should probably be easy, but I don't see yet how to prove it. Given is that $y=X_1\beta_1 + X_2 \beta_2 + \epsilon$. It holds that $X=\begin{pmatrix} X_1 & X_2 \end{pmatrix}$. Furthermore, it is given that $X_2'X_2=cX_1'X_1$. Now show using the formula $\operatorname{Var}(X'X)^{-1}=\sigma^2 V$, where $V=(1/(1+c))(X_1'X_1)^{-1}$. What I did so far: $X'X=\begin{pmatrix} X_1' \\ X_2'\end{pmatrix} \begin{pmatrix} X_1 & X_2\end{pmatrix}=\begin{pmatrix} X_1'X_1 & X_1'X_2 \\ X_2'X_1 & cX_1'X_1\end{pmatrix}$. However, I think it would be hard to invert this matrix. Any help is appreciated.


title: Type of ANCOVA design in R

text: I have the dataset that contains one continuous measure of plant leafs (PC1 scores of individual leafs), that is dependent variable and three others that I need to use as independant, explanatory. One is every leaf area, other is the individual and the third is the leaf number in sequence for individual. As every individual has certain and unequal amount of leafs I can't use simple ANCOVA. I guess there is a way of nesting the effect of leafs number in the effect of individual but I don't know how. I want to answer the question whether individuals are affected by size the same way? my google drive with the entire dataset in pdf. It is noticeable that one individual can have variable number of leafs, but also that the sequence of leafs is regular for every individual, from the lowest to highest. First column is leaf size, other is individual and then the leaf class followed by 4PCs.


title: How to rank criteria based on survey response?

text: I have made a survey to weight the importance of different features of a credit card using a 5 point Likert scale with scores (very important = 5...neither important or unimportant = 1). I have 8 criterions, each having a similar question, such as the one below: How important is the cash back feature on a credit card?               Very Important     Important     Slightly important     Not very important     Neither important or unimportant      All 8 questions follow the same format and wording, each associating with a different credit card feature. After conducting 10 surveys, what must I do to then determine which feature is the most important, in order to compute weights for each feature that equals to 1. Bear in mind that I am not a mathematician or a statistician. Thanks!


title: Using predict_proba with sklearn's multiclass SVC

text: I'm using python's sklearn for multi-class classification (SVC) When using the predict method, i get very high scores with my dataset, However, I want to plot ROC curves for each of my classes. That is, I would like to reduce the problem to a in_class/out_of_class problem for each of the classes. For that i resorted to the predict_proba method of the SVC. However, I find no correlation between the probabilities given and the predictions. For instance, for a 5 class classification problem, I may get a prediction of the 1st class, but get a probability vector - [ 0.1, 0.2, 0.5, 0.3, 0.0 ]. The 1st class did not get the highest probability. Does anyone know how the SVM uses its decision function to make a prediction, or how the predict_proba works on a multi-class problem? Thanks


title: Difference-in-difference with interaction terms

text: I am looking for an empirical paper that perform a difference-in-difference regression where the whole did regression (treated+post+post*treated) is interacted with a continuous variable (or a categorical one)... Do you have an idea of such a paper ?


title: Fitting a Non-Central t-Distribution with Location and Scale Transformations

text: I am trying to fit a distribution function to empirical observations that have the following properties:   1. Non-zero mean   2. Non-unit variance   3. Heavy tails   4. Asymmetric about the mode I am considering performing this task by fitting a non-central, generalized t-distribution, but I am struggling to find any literature online regarding estimation of the 4 parameters (i.e., degrees of freedom, location, scale, asymmetry parameter) from empirical data. Do closed-form estimators (MLE, method of moments, or anything else!) exist? Or is there perhaps a better distribution for this task, where such estimators do exist? I have looked at a 4-parameter logistic distribution - that looks like it would also do the trick, but I'm having similar problems finding the estimators that I require for this distribution. Many thanks!


title: Constructing an orthogonal polynomial regression model

text: My current model is: > $$y_i = a_0 +a_1x_i + a_2x_i^2 + \epsilon_i $$ With $E[\epsilon_i] = 0, Var[\epsilon_i] = \sigma^2$. I have the information: > $$\begin{matrix} i & \mathrm{Response} \, y_i & \mathrm{Covariate} \, x_i \\ > 1 & -1 & -0.5 \\ 2 & -0.25 &-0.25 \\ 3 & 1 & 0 \\ 4 & 1.2 & 0.25 \\ 5 & 2.6 > & 0.5 \end{matrix}$$ I need to propose an orthogonal polynomial regression model of the same order as my current one. So I thought that if I let $\zeta_1$ be the linear terms and $\zeta_2$ be the quadratic terms, then I want a model where $\sum \zeta_1 = \sum \zeta_2 = 0$. So I want two general forms of $\zeta$. So let us let: $$\zeta_1 = a + bx_i$$ $$\zeta_2 = a + bx_i + cx_i^2$$ Now I want to solve for $a,b,c$. But I'm stuck on how to solve for them. I can't just use normal simultaneous equations from here can I because everything cancels out leaving me with everything equal to $0$. What do I do?


title: Relative importance of a set of predictors in a random forests classification in R

text: I'd like to determine the relative importance of sets of variables toward a `randomForest` classification model in R. The `importance` function provides the `MeanDecreaseGini` metric for each individual predictor--is it as simple as summing this across each predictor in a set? For example:               # Assumes df has variables a1, a2, b1, b2, and outcome     rf <- randomForest(outcome ~ ., data=df)     importance(rf)     # To determine whether the "a" predictors are more important than the "b"s,     # can I sum the MeanDecreaseGini for a1 and a2 and compare to that of b1+b2?


title: statistical analysis of events in a given time interval

text: I am attempting to analyze biological data, to see whether the number of events in a given time interval is more/less than expected based on the overall frequency. How would one approach this? An example of how I would frame this: Out of 100 ms, 16/44 events occur in 15 ms, and 28/44 events occur in the remaining 85 ms. Do more events occur in the 15 ms interval than expected based on the overall frequency? And I guess the null hypothesis is that there is a random distribution of the 44 events over the 100 ms.


title: Designing Frequentist A/B Experiments

text: What parameters should one **pre-specify** before designing/running an experiment? (e.g. a two-sample, one-side, t-test) My understanding is that on establishes:   1. Null hypothesis   2. **Maximum** Significance level   3. **Minimum** Power   4. **Minimum** Effect size Is the choice of **maximum/minimum** in the lines above correct? Is it correct to say that one needs to establish all of the above to then estimate the **required** sample size? Or are there any other free parameters that one needs to establish first?


title: How to deal with non-normal heterocedastic data from a factorial experiment?

text: I ran an experiment with two factors each with two levels, 5 replicates each combination and one response variable. My data are non-normal and heterocedastic. Transformations didn't help. I ran a permutation ANOVA (`aovp - package lmperm - R`), but I am not sure this is the right way to go and I don't know what to do about a post hoc test, because I have a significant interaction. This is what my data look like:               F1 F2 BIOMASS     A C 0.66     A C 0.31     A C 0.88     A C 0.55     A C 0.81     A D 0.39     A D 0.21     A D 0.17     A D 0.15     A D 0.18     B C 0.00     B C 0.18     B C 0.05     B C 0.00     B C 0.02     B D 0.16     B D 0.04     B D 0.21     B D 0.15     B D 0.06


title: Can you use a spline function of the spatial coordinates to control for spatial autocorrelation?

text: > **Possible Duplicate:**   >  Why does including latitude and longitude in a GAM account for spatial > autocorrelation? I am interested in the effect of a predictor vector $X_i$ on a binary outcome $Y_i$, with corresponding spatial location $s_i$, so I've considered the spatial random effects model, \begin{equation} logit(P(Y_i=1)) = X_i \beta + Z(s_i) \end{equation} where $Z(s_i)$ is a stationary mean zero random process whose covariance structure is known up to a number of parameters $\theta$. It appears that this model is usually fit using Bayesian sampling algorithms, which can take a very long time to run when $n$ is large, since you have to sample from the joint posterior distribution of $\beta, \theta$ and the random effects. To circumvent the computational complexity and the expert tuning required for the MCMC, I've instead fit the model \begin{equation} logit(P(Y_i=1)) = X_i \beta + f(s_i) \end{equation} where $f(s_i)$ is estimated by a 2-D spline, which can be fit with the gam package in R in a matter of seconds. I have found in simulation studies that this model is effective at removing the spatial autocorrelation from the residuals and gives reasonable estimates for $\beta$. This leaves me wondering why this approach is not prominent in the literature. Why do users often opt for the computationally intensive sampling algorithms? Is there some drawback or something systematically lacking in this approach?


title: newbie question - regression vs composite normalised scores

text: Lets say I need to rank the "strength" of a bunch of objects based on a number of factors, as far as I can tell, there are two ways people seem to commonly do this : (a) calculate normalised scores (i.e. factor-avg(factor)/stddev(factor)), then sum the individual normalised results for each factor and rank those (b) do some sort of single or multivariate regression and rank based on that. How would you approach this sort of thing ?


title: Estimate of parameter $\beta$

text: Suppose I've a set of information, > >     id t  s   x          H0         h0        H0i         Hi Di       Eu > Elu >      1 2  1  53  0.01067033 0.01067033 0.01371871 0.01371871  2 1.039715 > 0.029300175 >      2 4  0  60  0.02148908 0.01081875 0.02856059 0.04227930  1 1.019138 > 0.009121439 >      3 5  0  45  0.04344224 0.02195316 0.05377424 0.09605354  2 1.038006 > 0.027655283 >      4 5  0  50  0.04344224 0.02195316 0.05506429 0.15111783  2 1.036866 > 0.026556713 >      5 6  0  52  0.05472922 0.01128698 0.07003181 0.22114964  2 1.035420 > 0.025161271 >      6 7  1  51  0.07763398 0.02290476 0.09887098 0.32002062  2 1.033386 > 0.023194494 >   If I want to find the estimate of $\beta$ from this, > $$ Q(\beta)=\sum_{i=1}^{G}\sum_{j=1}^{n_{i}}[\delta_{ij}\{\ln > h_{o}(t_{ij})+\beta'x_{ij}+E(\log > U_{i})\}-H_{o}(t_{ij})e^{\beta'x_{ij}}E(U_{i})] $$ where, $t_{ij}$ is the > $j^\mathrm{th}$ time point of the $i\mathrm{th}$ individual, $\delta_{ij}$ > is the status, $x_{ij}$ is the vector of covariates. Let initial $\beta$ be 0.02. How can we find the estimate of $\beta$ using R code?


title: "forgetfulness" of the prior in the Bayesian setting?

text: It is well-known that as you have more evidence (say in the form of larger $n$ for $n$ i.i.d. examples), the Bayesian prior gets "forgotten", and most of the inference is impacted by the evidence (or the likelihood). It is easy to see it for various specific case (such as Bernoulli with Beta prior or other type of examples) - but is there a way to see it in the general case with $x_1,\ldots,x_n \sim p(x|\mu)$ and some prior $p(\mu)$? EDIT: I am guessing it cannot be shown in the general case for any prior (for example, a point-mass prior would keep the posterior a point-mass). But perhaps there are certain conditions under which a prior is forgotten. Here is the kind of "path" I am thinking about showing something like that: Assume the parameter space is $\Theta$, and let $p(\theta)$ and $q(\theta)$ be two priors which place non-zero probability mass on all of $\Theta$. So, the two posterior calculations for each prior amount to: $$p(\theta | x_1,\ldots,x_n) = \frac{\prod_i p(x_i | \theta) p(\theta)}{\int_{\theta} \prod_i p(x_i | \theta) p(\theta) d\theta}$$ and $$q(\theta | x_1,\ldots,x_n) = \frac{\prod_i p(x_i | \theta) q(\theta)}{\int_{\theta} \prod_i p(x_i | \theta) q(\theta) d\theta}$$ If you divide $p$ by $q$ (the posteriors), then you get: $$p(\theta | x_1,\ldots,x_n)/q(\theta | x_1,\ldots,x_n) = \frac{p(\theta)\int_{\theta} \prod_i p(x_i | \theta) q(\theta)d \theta}{q(\theta)\int_{\theta} \prod_i p(x_i | \theta) p(\theta)d \theta}$$ Now I would like to explore the above term as $n$ goes to $\infty$. Ideally it would go to $1$ for a certain $\theta$ that "makes sense" or some other nice behavior, but I can't figure out how to show anything there.


title: Increasing the sample size does not help the classification performance

text: I am training a SVM classifier based on a given document collections. I started from using 500 documents for training, then I add another 500 for training, and so on. In other words, I have three training sets, 500, 1000, 1500. And the smaller training set is a subset of the sequential larger set. I validate the model against the same test set.               SampleSize  Precision   Recall      Accuracy      AUC       500       79.62%          67.49%      77.65%    0.854       1000       82.49%         77.94%       82.67%   0.890       1500       81.64%         78.08%       82.28%   0.888      The performance gets the best when we use 1000 training set. Looks like the extra 500 training samples for constructing the 1500 training set in fact does harm to the model. How can I explain this observation.


title: Importance of the power of t-test in the case when $H_0$ is rejected

text: I know that $\beta$ (type II error) is the probability that the null hypothesis is false, but it is not rejected. The power of the test is defined $1−β$. I still thought that the power of the test is important only in the case when $H_0$ is not rejected. However my statistical software compute the power of the test also in the case when $H_0$ is rejected. In my opinion information about power of the test in the case when $H_0$ is rejected is useless. Is it true? Or is information about power in the case when $H_0$ is rejected usefull for some kind of statistical decision.


title: Luck or skill in 3-person game

text: Two friends and I have been playing a game that is a combination of skill and luck. (as most games are). We assume that if the game was all luck and/or we all had the same skill level eventually the win% of all of us would approach 33%. So what sample size would we need to show that one of us was "better" at the game? That is if we played 50 games, what % of games (or how far from 33%) would someone need to demonstrate skill vs. luck (perhaps at a 95% or 99% likelihood). What if we played 100 games or 200 games? I would assume we would need a smaller percentage difference from 33% as the sample size went up.


title: Summary statistic for 'average' conditional probabilities

text: I have the following scenario: Let $X_i$ denote some event. Now I have the following calculated probabilities. $P(X_1 \cap X_2 \cap \cdots \cap X_9 | X_{10})$ $P(X_1 \cap X_2 \cap \cdots \cap X_{8} \cap X_{10} | X_{9})$ $\vdots$ $P(X_2 \cap \cdots \cap X_{9} \cap X_{10} | X_{1})$ In words, the above expressions represent the (joint) probability of the 'remaining' events given that a particular event has occurred. To calculate these probabilities, I know the underlying probability distribution that describes the entire system, i.e., $p(x_1, x_2, \cdots. x_{10})$, so for example, calculating the probability $P(X_1 \cap X_2 \cap \cdots \cap X_9 | X_{10})$ just simply requires me to find $\frac{P(X_1 \cap \cdots \cap X_{10})}{P(X_{10})}$ where both the denominator and numerator can be calculated by integrating over certain regions of the probability density function $p(x_1, x_2, \cdots. x_{10})$. My question is, I wish to find one summary value that describes the 'average' probability of this system, e.g. say $P(X_1 \cap X_2 \cap \cdots \cap X_9 | X_{10})= 0.4$, $P(X_1 \cap X_2 \cap \cdots \cap X_{8} \cap X_{10} | X_{9})= 0.3$ etc, how can I "combine" these values $0.4$, $0.3$, etc into one value that describes the "average" probability of this system? My initial method is just to take the arithmetic average of each value, but that is mathematically incorrect since conditional probabilities aren't summable (except when they are conditioned on the same event). So are there any other measures/techniques I can use to somehow "combine" these single probabilities into "one" value? I.e., some kind of summary statistic that describes the average probability?


title: Separating mixture with little a priori knowledge

text: I have a dataset (time series data) of measured signal power values from a radio receiver. The data does not originate from a controlled experiment. I have limited knowledge of the underlying processes and distributions. My observations are the result of one or more processes. In any case, there is a background noise process. Additionally, there can be one or more processes belonging to signals of radio transmitters. Here are three examples (histograms). In the first case, there is only noise. The second and third case, there is a real mixture of noise and more than one radio transmitter. ![only noise example](http://i.stack.imgur.com/7Pj0Z.png) ![mixed example 1](http://i.stack.imgur.com/GV7S8.png) ![mixed example 2](http://i.stack.imgur.com/Hooxs.png) My problem is the following. I want to estimate the fraction of samples that belong to the background noise process. There can be zero to many "non-noise" processes and I have no a priori information about this number. In general, the sub-populations do not follow the normal distribution (the x-axis in the graphs above is actually logarithmic, note the "dB".) I have no a priori information about their mean and variance. The dataset is huge. Therefore, I am looking for an approach that does not need manual adjustment. So far, I worked with a simple threshold to separate noise from the rest, but this very inaccurate, especially when there is a signal process with mean very close to the noise. Can you point me to methods which can be applied here? **EDIT** : Further reading makes me think that I may want to:   * Derive a mathematical description of the _background noise_ distribution (curve) from my empirical data.   * Fit this curve to the sample data.   * Calculate the ratio of the area under the fitted curve and the total number of samples. The question is then how to fit the curve? The algorithms I found so far either need the number of components a priori and/or have limitations about the distributions they assume.


title: Two-way MANOVA vs. two-way ANOVA on a composite variable

text: Is using a 2 way MANOVA to check whether height and weight differ across age groups and gender, the same as using a 2 way ANOVA to check whether BMI (wt in kg/ht in m$^2$) differ across age groups and gender?


title: Divergence measure of two classifiers' performance?

text: I have two classifiers built with the same data. How can I measure divergence of these models? I found something like DIC but I don't know how to calculate this in R?


title: Using Eta appropriately

text: I have a nominal variable and a continuous (interval) variable. SPSS crosstab includes an option for Eta. Is it appropriate to use eta to test effect of the binary nominal variable on the continuous variable? How would I interpret the result such as Eta = 0.04? Thanks.


title: Using Normal Q-Q Plot to detect outliers in non-normally distributed values

text: I have values that are not normally distributed since most (132) of the (140) values are 0 or 1. There are a few that are <10, one about 80 and one about 650. Is it valid to use a normal Q-Q plot to identify the outliers?


title: calculation of composite scores from factor analysis

text: We're designing a survey where we ask students to indicate their level of agreement (likert scale 1 - Strongly agree....5 - Strongly disagree. We ran factor analysis on the data and obtained 4 factors. In order to get composite scores for each factor, this is what we did: For example, If factor 1 was made up of 3 items, each with raw scores of 3,4 and 5, we took the mean of 3 scores and then a z-score of the mean to get standardized scores for further analyses. Now if we want to divide the construct being measured into high and low points, would z-scores > 0 automatically imply a high measure of the construct? In the original scale, if 3 was the 0 point, values above 3 would indicate disagreement (low measure of the construct). So does the same thing translate to the z-score as well? Would z-scores above 0 indicate low measures? Thanks


title: How to interpret the results from random coefficient Poisson data analysis

text: I have fitted random coefficient Poisson analysis in R. I have obtained the following results: Generalized linear mixed model fit by maximum likelihood ['glmerMod'] Family: poisson ( log ) Formula: frequency ~ 1 + cc + ageveh + make + (1 | AREA) Data: x                 AIC       BIC    logLik  deviance       1359.1477 1389.7370 -672.5739 1345.1477 Random effects: Groups Name Variance Std.Dev. AREA (Intercept) 1.323 1.15 Number of obs: 584, groups: AREA, 8 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.12902 0.44432 -0.290 0.7715 ccL 0.05656 0.12371 0.457 0.6475 agevehO 0.02136 0.09264 0.231 0.8177 make2 -0.45454 0.20632 -2.203 0.0276 * make3 -0.31799 0.21422 -1.484 0.1377 ## make4 -0.29708 0.14469 -2.053 0.0401 * Signif. codes: 0 ‘ ** _’ 0.001 ‘** ’ 0.01 ‘_’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects:                   (Intr) ccL    agevhO make2  make3      ccL 0.052 agevehO -0.179 -0.232 make2 -0.171 -0.007 -0.001 make3 -0.156 0.022 -0.078 0.366 make4 -0.300 -0.235 0.167 0.544 0.522 However I am unable to interpret the results.


title: Algorithm for online handwriting recognition

text: Is there any specific algorithm for online handwriting recognition? The algorithm should recognize non-cursive and cursive handwriting. I know there is already a similar post on stackoverflow.com, but it was posted 5 years ago, and I did not find an answer that satisfies my need. I assume within these five years, advanced algorithms for online handwriting recognition could have been developed. Any help would be greatly appreciated.


title: A seeming paradox with rational agents not coming to the same conclusion given the same data

text: So one day after a tasty dinner full of bananas, an idea comes to your mind (you are person A) - "What if eating bananas can cure cancer?". Being a scientist at heart, you conduct a double blind study and your data shows you that bananas cure cancer with probability at least 99.99%. "WOW!" - you think, - "I am onto something important here!". You pack your stuff and hastily head towards the central organization for cancer research in your country. "People!", you exclaim entering the building, "I found the cure for cancer!". They direct you to the head of the organization (person B). You show him your data, and he agrees that your conclusion based on your data is correct. "But, you see," - says the guy, - "this year we got at least 100,000 people like you coming to us and reporting their findings of curing cancer with various kinds of fruits. And we expect at least some of them to have data showing that they cure cancer with probability of 99.99% purely by chance. Therefore I am unconvinced, and you should conduct further research.". "So what?" - you say, - "I didn't conduct all those experiments, I only conducted one experiment, and it showed me that bananas cure cancer with probability 99.99%. Therefore I am bound to believe that bananas cure cancer with probability 99.99%. All the other studies used different fruits, and none of them used bananas, so they are irrelevant for this matter". Assume that both person A and person B are absolutely rational and they completely trust each other. It is generally assumed that two rational agents given the same knowledge must come to the same conclusions. In this case though it appears to me that there is nothing that person B can tell person A that would lead person A to accept the person B's position. What will happen I think, is that person A will accept that person B should indeed hold the position that it is not clear whether bananas cure cancer or not. And person B will accept that person A is justified in believing that bananas cure cancer with probability 99.99%. So this appears to be a paradox to me: from one side rational agents given same information must come to the same conclusion, and from the other hand we have here that rational agents have the same information but come to different conclusions. Do you see any resolution to this?


title: Sampling from a fixed population

text: Here's a real basic question. I'm trying to teach myself a bit of stats with Verzani's Using R for Introductory Statistics. In question 5.13 he asks: A sample of 100 people is drawn from a population of 600,000. If it is known that 40% of the population has a specific attribute, what is the probability that 35 or fewer in the sample have that attribute. Now, I guess you're supposed to reason that the population is sufficiently large that assuming independent Bernoulli trials is close enough. Then, you get your answer like this:               > pbinom(35,100,0.4)          [1] 0.1794694      My question is this. How would you go about answering a question like that without assuming independence, say if the population was smaller. I'm sure it'll become obvious after I read more. Just trying to make sure I'm not missing something. Sorry for the intro level question. Thanks!


title: Bayesian user survey with a credible interval

text: I work on a couple of library websites and we want to know who are users are so we can make decisions about how to design the site and focus efforts. These are the classes of users we are interested in: "Please select the answer that best describes yourself:"   * k12-teacher-librarian   * k12-student   * college-student   * grad-student   * faculty-researcher   * genealogist-researcher   * archivist/librarian   * other [user can specify or leave blank] We want to know the proportion of our users who fit into each category, and we want to have some sense of how accurate our numbers are. We also want to know if/how these proportions change during different times of year. My first question is "is this a problem where I can use Bayesian stats?" We ran a pop-up survey that went on every page of both our websites for two weeks about 6 months ago. It used a flash cookie to mark when the survey had been taken so the users would only get this once. This got logged into google analytics so we could see the % of use of the sites from each type of user. However, our websites are used a lot in library computer labs where they are set up with flash on a read only filesystem which ended up with the survey showing up on every single page since we could not really write to the flash shared object. Can I use the data from the "shotgun" pop-up survey as my priors; and then prospectively run the pop-up survey for one out of every 1,000 unique visitors to each site? I'd have 6 to 8 new survey results coming in a day that I'd like to use to continually update a Bayesian model of who our users are. The thing is, I have not had a stats class in 20 years and then it was classical stats. I've been googling around about discrete random variables and reading http://www.greenteapress.com/thinkbayes/ but I really have no idea what I'm talking about. Am I on the right track with this idea (if so, any pointers how to proceed? I do work in a library so I can look up any reference you can point me at) or does this sound like hogwash/gibberish? * * * update: Here is the data from the survey we ran in December 2011               Website A:      65    k-12 teacher or librarian      71    k-12 student     532    college or graduate student (we want to split this going forward)     307    faculty or academic researcher     369    archivist or librarian     234    genealogist or family researcher     584    other          Website B:     133    k-12 teacher or librarian     280    k-12 student     445    college or graduate student      82    faculty or academic researcher      55    archivist or librarian      68    genealogist or family researcher     220    other


title: A re-formalization of a conjugate prior?

text: It is quite easy to prove that if $p(\theta)$ is a conjugate prior to some likelihood then the following: $$q(\theta') \propto p(\theta)I(\theta \in A)$$ where $A$ is a subset of the parameter space and $I(r)$ is 1 if $r$ is true and 0 otherwise -- is also a conjugate prior. Basically, if we take a conjugate prior and re-normalize it to a subset of the parameter space, we still get a conjugate prior family. Is this mentioned somewhere that I can refer to?


title: How to control for common method bias?

text: I am doing CFA to test for common method bias using an unmeasured latent construct method. The average common variance I am getting from the common factor is about 31% which is high (as should be less than 25%). I checked for differences in p-values of the trait loadings (between CFA model with ULC and without ULC) and p-values become insignificant for only one scale. What can I do to control the effect of common method variance statistically in this case? I am doing a moderation analysis and I have 6 variables


title: Random sample of population sub sample size

text: I am generating models and each model is a random sample of the total model population. It is recommended that I generate 30,000 models and cluster taking the top 5 to 10 clusters to reach the native state. I would like to adjust the 30,000 number lower but still remain within the probable range of the native structure. If $\Pr(\text{Native}\mid 30000)=1$; what is $\Pr(\text{Native}\mid X)$? I read about Bayesian statistics but I don't have all the pieces to put together Bayes' theorem for this problem. I also read about frequentist approach but I haven't used it before. "Cluster analysis or clustering is the task of assigning a set of objects into groups (called clusters) so that the objects in the same cluster are more similar (in some sense or another) to each other than to those in other clusters."


title: Fixed-sum question

text: I have a set of 10 items each with a weight between 0-1. A person can rate each of these 10 items with a score between 0-100, but is given, say, only 500 points to allocate among these 10 items The person does not know the weight of each item and the total score will be the weighted average. A person could rate each item as 50 and another person could rate the top 5 items with the highest weights with all 100's. The weights are such that the difference in total scores received by rating all items with 50, and only the top-5 highest weighted items with 100s and the remaining 5 items with 0, is not significant given that the spread in the weights (all between 0-1) is not high. Given that a person can randomly rate all the items with 50's and still receive a total score very close to the highest possible score obtained by rating the top-5 weighted items with all 100's, is problematic. Are there any existing alternative approaches/methods/algorithms to solve such a problem so that a person who rates all items 50 will get a score far away from the maximum possible score?


title: Estimating one logit vs. two logits using 2 dependent variables (i.e. restricting the two logits to have the same coefficients

text: I have a data set where a person has no disability, a mild disability, or a severe disability. I want to determine of the model for mild and severe disabilities are the same. I estimated three logits using the same set of independent variables (X's) In the first model the dependent variable Y1 = 1 if they have any disability (mild or severe) and 0 otherwise Then I estimated two other models. The first has the dependent variable Y2=1 if they have a mild disability or 0 otherwise. And then I estimated a logit with Y3=1 if they have a severe disability and 0 otherwise. (In other words Y1=Y2+Y3) How do I use the the results to test if anything is gained by estimating the two models separately, or if it is better to estimate the first model, thus restricting the coefficients on for models 1 and 2 to be the same


title: Winning percentage -- logistic regression or linear regression?

text: Suppose I want to know whether a baseball team's winning percentage the previous season or the Pythagorean expectation of that team from the previous season is a better predictor of next season's win percentage. (I take percentages because the number of games played per season has changed in the past, so it seems easier). For each of these independent variables, there's good reason to suspect their will be a linear relationship with the dependent variable, so linear regression seems like it would work. The dependent variable is always within the [0,1] interval, so logistic regression also seems like it would work. But the dependent variable isn't exactly a probability (well, I suppose it's the probability that the team beats an unknown opponent team), and it also never is actually very close to 0 or 1 (it's pretty much without exception between .3 and .7). So, with all this in mind, which would be the more natural method to use, linear regression or logistic regression? Are they both valid approaches?


title: Prediction and Tolerance Intervals

text: I have a couple of questions for prediction and tolerance intervals. Let's agree on the definition of the tolerance intervals first: We are given a confidence level, say 90%, the percentage of the population to capture, say 99%, and a sample size, say 20. The probability distribution is known, say normal for convenience. Now, given the above three numbers (90%, 99% and 20) and the fact that the underlying distribution is normal, we can compute the tolerance number $k$. Given a sample $(x_1,x_2,\ldots,x_{20})$ with mean $\bar{x}$ and standard deviation $s$, the tolerance interval is $\bar{x}\pm ks$. If this tolerance interval captures 99% of the population, then the sample $(x_1,x_2,\ldots,x_{20})$ is called a _success_ and the requirement is that 90% of the samples are _successes_. Comment: 90% is the _a priori_ probability for a sample to be a success. 99% is the conditional probability that a future observation will be in the tolerance interval, given that the sample is a success. My questions: Can we see prediction intervals as tolerance intervals? Looking on the web I got conflicting answers on this, not to mention that nobody really defined the prediction intervals carefully. So, if you have a precise definition of the prediction interval (or a reference), I would appreciate it. What I understood is that a 99% prediction interval for instance, does not capture 99% of all future values for _all_ samples. This would be the same as a tolerance interval that captures 99% of the population with 100% probability. In the definitions I found for a 90% prediction interval, 90% is the _a priori_ probability given a sample, say $(x_1,x_2,\ldots,x_{20})$ (size is fixed) **and** a single future observation $y$, that $y$ will be in the prediction interval. So, it seems that both the sample and the future value are both given at the same time, in contrast to the tolerance interval, where the sample is given and with a certain probability it is a _success_ , and under the condition that the sample is a _success_ , a future value is given and with a certain probability falls into the tolerance interval. I am not sure if the above definition of the prediction interval is right or not, but it seems counterintuitive (at least). Any help?


title: Total probability of any event

text: We have an application where we're trying to quantify the total risk(hitting traffic, say) of a trip in a car where probability is a function of time and a series of algorithms: ![Our probability graph with error bars](http://i.stack.imgur.com/nMP78.png) How can we quantify the 'total' risk of the trip while respecting the amount of time the car is exposed to higher risk? My gut says it's a matter of an integral but my stats isn't strong in the least and think I might be missing something obvious. The y-axis is the probability that an event occurs and the x-axis is the time elapsed, in minutes. The gray area is irrelevant.


title: Interpreting the results of a binary logistic regression

text: This is a basic question. I have been handed a binary logistic regression. The model has significant terms, but the goodness of fit tests indicates the logit model is not appropriate. The author of the study indicates that the goodness of fit data does not invalidate the relationship between the dependent variable and the predictors, only the ability of the model to accurately predict outcomes. the argument is that since we were only interested in verifying a relationship and not the magnitude, the result is conclusive. I'm skeptical of this. Wouldn't it be more appropriate to say that the lack of fit does not necessarily invalidate the relationships? With a different link function, couldn't the observed and expected counts change enough to move some insignificant term to significance or vice versa?


title: Imputation with Random Forests

text: I have two questions on using random forest (specifically randomForest in R) for missing value imputation (in the predictor space). 1) How does the imputation algorithm work - specifically how and why is the class label required for imputation? is the proximity matrix which serves to weight the average value to impute a missing value defined separately by class? 2) If the class label is needed to impute missing values - how can this be used to impute missing values for new data that you are trying to predict?


title: Cluster analysis on panel data

text: I have a panel data set (country and year) on which I would like to run a cluster analysis by country. My data set has around 20 variables. Here's a summary for my panel data: `panel variable: country (strongly balanced) time variable: year, 2010 to 2013` Running a kmeans cluster analysis on 2013 data only is pretty straightforward. But how would you do the analysis considering all observations in the 2010-2013 period? Is k-means clustering an appropriate approach? Here's what I ran in Stata for 2013 only: `cluster kmeans var1 var2 var3 var4 var5 var6 if year==2013, k(4) name(test1)` Thanks!


title: How to choose a sample size for two group repeated measures anova design given power, alpha, and effect size?

text: How would I go about computing the sample size for repeated measures ANOVA design, given power, significance and the required effect size (2 groups, I know the mean difference that I want to be able to catch and approximate standard deviation)?


title: What do you do after testing if two non-normal distributions differ?

text: This is a follow up question to "How do I test if two non-normal distributions differ? I have 13 distributions, ![X vs Y](http://i.stack.imgur.com/eLKCE.png). According to the previous question, I should do the following.   1. Kruskal-Wallis test with 13 samples.   2. Mann-Whitney U pairwise.  I found that all p-values are practically 0. At first, this did not seem right. For example, distributions A and B look very similar. (Subquestion: Should I use `alternative = "greater"` in R for example? I would like to roughly argue that D is greater than A.) Now, based on another question, I believe that this is due to the fact that my sample sizes are so large. Based on that same question, it appears that I need to calculate the effect size. I followed these instructions. I obtained tiny effect sizes, for example, 0.00005. How can I interpret these values? I obtain "tiny" effect sizes even when comparing A and D, which to me, look like very different distributions. To address `@whuber`'s comment, I would like to state that D > L > K > J > ... > [A,B,C]. In other words, that those in D are generally larger than those in L. In the case of A,B,C, I would have to say that there is no significant difference. But again, I would like to an argument based on statistics not just visual observation.


title: Does differential geometry have anything to do with statistics?

text: I am doing master in statistics and I am adviced to learn differential geometry. I would be happier to hear about statistical applications for differential geometry since this would make me motivated. Does anyone happen to know applications for differential geometry in statistics, please? Thank you!


title: Visualizing counts that have some outliers

text: I am running a project which compares citation counts across datasets per year, and then wants to plot the general change in the frequency in which a given source is cited. Here are two example datasets: 1999 - SOURCE A               {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 23, 70, 0, 0, 0, 0, 0, 0, 0, 32, 0,      0, 0, 25, 0, 3, 13, 0, 1, 0, 0, 0, 0, 11, 0, 0, 0, 0, 0, 0, 0, 0, 8,      0, 0, 6, 0, 0, 0, 1, 0, 22, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 2, 0,      21, 6, 10, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0}      2000 - SOURCE A               {0, 0, 0, 0, 0, 4, 400, 82, 0, 0, 2, 0, 0, 0, 2, 6, 0, 0, 12, 24, 0,      0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 0, 142, 0, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 22, 0, 0, 0, 0, 0, 10, 0, 4, 0, 0, 0}      In this case, when I plot the difference between 1999 and 2000 in a simple line-chart, 2000 seems so much bigger thanks to the data points that have 400, 142. Is there a good way to visualize this kind of data that stays true to outliers, but still shows the general trend that's going on (if any)? Thanks in advance.


title: How to detect step changes in GPS time-series data?

text: The graph below shows GPS heading data (sampled every second) and I am trying to find the best way to detect (right/left) turns in the data. Appreciate suggestions for algorithms/methods for it (perhaps step detection)? ![GPS heading sampled 1 sec apart](http://i.stack.imgur.com/nFMPt.png) While turns typically result in sharp change in heading value, heading might also change gradually due to road curvature which should not be detected.


title: Aggregating or Pooling Multinomial Logistic Model Coefficients

text: I'm trying to recreate the MNL analysis like this paper or this one. In Excel I have a sheet with my choice column which contains the market share on each itinerary. Alternatively I also have the number of passengers on that itinerary that could be used instead as the choice variable. I then, have several variables that describe the utility of the given itinerary that are similar to the ones in the papers above as my independent variables. The problem I'm having is that rather than spiting out the coefficients for my independent variables for the entire dataset like shown on page 22 of the first paper or page 18 of the second paper, I'm getting the variables for each choice variable in a huge list. Currently I'm trying to use the `mlogit` package in `R`. I have also tried `XLStats` with similar challenges. Accord.NET seems to produce the output I like but it appears to be doing some kind of binomial logistic regression and not MNL. **simplified data example**               Response variable: Market Share            Predictor 1: Origin/Dest/Carrier (nor used in model)     Predictor 2: IsNonstop?(binary variable)     Predictor 3: NumberOfSeats      * * *               50% DEN/JFK/AA 1 5000      20% DEN/JFK/UA 0 3000      30% DEN/JFK/DL 0 8000      25% DFWLAX/F9 0 10000      75% DFW/LAX/WN 1 25000       **desired output**               IsNonstop? 2.30      NumberOfSeats .12       **actual mlogit output**               50% IsNonstop 1.9 NumberOfSeats .45      20% IsNonstop 1.4 NumberOfSeats .015     30% IsNonstop -41.9 NumberOfSeats .85       etc....


title: Issue with controlling confound in multinomial regression analysis; different results when removing kids on meds

text: I examined the influence of ADHD on abnormal bodyweight in a very large, national sample of children. In my multinomial regressions, I controlled for several specific confounds, which have been shown to be associated with both ADHD and bodyweight in previous studies. One of these confounds was ADHD medication. Although it was significantly associated with underweight (and not obesity; both relative to normal weight), the results I got were different when I reran my analyses after excluding children on both ADHD medications and "other" medications. When I was statistically controlling for medication use, ADHD meds increase odds of underweight, and "other" meds decreased odds of underweight (neither impacted obesity), and the relationship between ADHD and underweight was significant in the model. However, when I removed kids on these meds, the association of underweight and ADHD became non-significant (very, not even "marginally" for those who like shades of grey). Are there particular problems that can arise when trying to control for other variables in multinomial regression, and which may have led to this situation? Clearly, the model "controlled for" medication use, since both types of medications were significantly associated with the outcome (odds of being underweight vs. normal weight), but it seems like it did not completely account for the impact of medication use in the relationship between ADHD and abnormal bodyweight.


title: Why does the log likelihood need to go to minus infinity when the parameter approaches the boundary of the parameter space?

text: In a recent lecture I was told that, in order for the maximum likelihood estimate to be valid, the log likelihood needs to go to minus infinity as the parameter goes to the boundary of the parameter space. But I don't understand why this is essential. Suppose the log likelihood goes to some kind of asymptote. Then the parameter that maximizes the likelihood is still the maximum likelihood estimate, right?


title: Cointegration but no Granger-Causality found

text: I have found Cointegration based on Engle/ Granger and Johansen. However, Granger-causality is rejected for both variables. How is that possible? According to theory, if x and y are I(1) and cointegrated, x is Granger causal to y and/or y is Granger causal to x. However, Granger-causality has been rejected in my bivariate case, despite their cointegration relationship. Did I understand it correctly that there has to be at least one granger causality flow in a bivariate cointegrated system? Thank you for your answer! Applying an VECM, I get the following results: with only the -0.022460 being significant...                   Vector Error Correction Estimates            Date: 10/28/13   Time: 23:58             Included observations: 1113 after adjustments            Standard errors in ( ) & t-statistics in [ ]                  Cointegrating Eq:   CointEq1              CAD(-1)             1.000000              NATGAS(-1)          0.067366                            (0.02646)                            [ 2.54615]               C                  -0.077093               Error Correction:   D(CAD)  D(NATGAS)          CointEq1        -0.022460   -0.006601      (0.00514)   (0.01384)     [-4.37213]  [-0.47714]          D(CAD(-1))  -0.054710    0.029241       (0.02998)  (0.08073)     [-1.82508]  [ 0.36220]          D(CAD(-2))   0.035656    0.101838      (0.02996)   (0.08070)     [ 1.18998]  [ 1.26200]          D(NATGAS(-1))   -0.004642   -0.077700      (0.01120)   (0.03016)     [-0.41449]  [-2.57591]          D(NATGAS(-2))    0.004712    0.056858      (0.01120)   (0.03016)     [ 0.42067]  [ 1.88491]          C    0.000176   -0.000850      (0.00019)   (0.00051)     [ 0.92332]  [-1.65571]           R-squared   0.022437    0.011948


title: Is there anything special about Gamma distribution with the shape parameter k=e?

text: Is there any unique property of $\mathrm{Gamma}(k=e, \text{ scale})$ or a Negative binomial distribution with $r=e$? Here, $e$ is Euler's number, $e \approx 2.71828$. The reason I'm asking is that one of the variables in my computer simulations can be fitted by $\mathrm{Gamma}(k=e, \text{ scale})$ or by $\mathrm{NB}(r=e,\ p)$ very robustly. That makes me wonder if there is something special about this process and this particular value of the shape parameter might hint to that.


title: Do Bayesians ever argue there are cases in which their approach generalizes/overlaps with the frequentist approach?

text: Do Bayesians ever argue that their approach generalizes the frequentist approach, because one can use non-informative priors and therefore, can recover a typical frequentist model structure? Can anyone refer me to a place where I can read about this argument, if it is indeed used? EDIT: This question is maybe phrased not exactly the way I meant to phrase it. The question is: "is there any reference to discussion of the cases in which the Bayesian approach and the frequentist approach overlap/intersect/have something in common through the use of a certain prior?" One example would be using the improper prior $p(\theta) = 1$, but I am pretty sure this is just the tip of the tip of the iceberg.


title: Effect of mtry and trainingset size in Random Forest

text: I've plotted the learning curves below, using different RF trainers. The training set is very small and with few features (it's the popular _Titanic_ dataset). At pre-processing stage, I just created a set of dummy vars from my original factor features. I then used the AUROC to choose the best value for mtry. But I'm not sure about the best interpretation of these learning curves. ![RF learning curves](http://i.stack.imgur.com/pTB84.png) The left one shows the missclassification error for varying training set sizes. The training set accuracy decreases while increasing the training set size, which is something I really expected. For the same reason, I expect less test error while increasing the training set size, which is what I got (more or less). The right plot shows the missclassification error when varying the parameter 'mtry'. Here, I got something really unexpected. I tuned the model to mtry=4 using CV. While I was expecting a decreasing missclassification rate in training set while increasing mtry (probably due to overfitting), I expected a very low value (if not a global minimum) on mtry=4. Even more weird, on the test set: it seems to be fairly constant/slowly decreasing, with a little spike (noise?) on mtry=4! This is really an unexpected result! What is causing me to observe these trends?


title: Fitting NLS algorithm

text: I have the list of time series. I am fitting these series with the formula $$y=ax^2\exp(b*x)$$. It should be noted that parameter $b$ in the formula must be negative as this reflects the behaviour of time series. An example plot for the fitted regression output using `nls()` in R is: ![Fitting with nls\(\)](http://i.stack.imgur.com/OoMEZ.png) Initially, I was faced with the problem of choosing start values for $a$ and $b$. I chose such values so that the fitted value evaluated at the maximum for $y$ is the same. However with some time series I faced with the next problem. Let's consider this time series: ![Increasing time series](http://i.stack.imgur.com/7VfnM.png) When applying my methodology, `nls()` gave me a positive value of $b$ which is unsatisfactory for me although the fitted values were highly concordant. However when I tried to introduce boundaries in the `nls()` function such as `upper(a=10,b=-0.0001)`, I found that no matter what value of $b$ I put, the resulting value of $b$ does not converge to any stable value and the fitted curve is highly discordant. Numerically, the parameters were consistent with what I expected to find, but the fitted values were uncalibrated. Can anybody give any suggestions of how can I deal with such time series without changing the general formula? Thanks!


title: How can I estimate optimal cut-off using linear regression-models

text: To build two linear regression model, (dependant var : B, independant vars : A1, A2, A3) I have to set the cut point of A1. (high A1 and low A2) I want to pick up the model*s*(a model for high A1, and the other for low A2) with with the lowest residual sums of square*s* In summary, How I can get the optimal cut-off point to get best linear regression models for each groups? The article I am trying to replicate is below, that description is all. A is dep Var, and B is ind var. "Considered two linear regression models, one for subjects below a certain level of A and the other for subjects above that level. To determine the specific cutoffs, we fitted the two linear regression models described above and calculated the sums of squares of residuals (=observed B - estimated B) from the two models for each level of A. The models with the lowest residual sums of squares were our best models, and the corresponding level of A were defined as the optimal cutoff values." Url for the article is : http://link.springer.com/article/10.1007/s00223-012-9669-3


title: Probability to get equal cars

text: A famous car production company produces 3 different models of cars, available in 5 colors and 3 possible CC(1900 cc,2000 cc and 2200 cc). suppose that every combination of model, color and CC has equal probability to be selected by a client:   * If in one week 6 cars are sold, what is the probability that at least two cars are equal(same color,cc, and model) My answer: Total possible combinations are 45^6, instead of finding "at least" two cars i would find 1-(Probability to get zero equal cars + probability to get one equal car) :               P(get 0 equal cars)= 44^6     P(Get 1 equal cars) 44^5      So the answer is: $1-(44^6+44^5/45^6)= 1-0.893= 0.107$ Is my answer correct? thanks


title: 3 way validation - model results very different between validation columns

text: I have a dataset of approx 1000 members and 100 variables modelling is multivariable using mainly partition methods and also some neural net experimentation I have been experimenting with 3 way validation using both 60-20-20 and 70-15-15 the problem i have is this: Using 3 way validation and running a model generated using one specific validation set then against with differing 60-20-20 set produces quite different results, sometimes contributing predictor variables will change significantly and also r^2 for test portion can vary between 0-0.5 this says to me i do not have enough data entries in this dataset note that putting all 100ish predictors in results in only 10 or so being significant in any one model with the upper most usually contributing even between different validation sets (but not always in simlar arrangements of contribution) I think i am at the point that if i exclude some colinear variables i can viably use way as opposed to k-fold, [k-fold always splits too far requiring manual pruning]. am i indeed in the area of being able to use 3 way or not?


title: Algorithms/methods for detecting regime changes in data rates

text: I'm looking at multiple live data streams which I process on a daily basis and would like to monitor in real time. The way I did it at first was to introduce a somewhat arbitrary upper bound for each of the streams..the monitoring would look at time elapsed between consecutive stream items and if it exceeds its threshold it would alert on possibly missing data on that stream. The problem with this approach is that the streams are very active at some points of the day and rather quiet at other points, so a hard coded limit ends up being unnecessarily high to avoid false alarms during the tranquil periods. The second way I reworked it was to introduce bucketed data rate thresholds for each data stream. So, for example, something like: Stream A 4am threshold, Stream A 7am threshold, Stream A 5pm threshold...and the threshold used throughout the day is dependent on the time. This is more reasonable but the problem I have now is in coming up with appropiate times for each stream in a better way than arbitrarily choosing them for each data stream. For a given stream, I envision stepping through the data and generating "reasonable" buckets along with their threshold (for now I'm using twice the maximum time between items because their distribution is not approximately normal so I haven't looked at how to get a threshold time which would correspond to a certain low probability of occurrence) but I'm wondering where to start for this? And is this similar or identical to a well studied problem with solutions I might want to use?


title: Posterior distribution

text: Suppose $X1,..,X4$ be iid from pdf $f(x|\theta)=\frac{1}{\theta}$ ,for $0<x<\theta$. The prior distribution is $\pi(\theta)=\frac{2}{\theta^3}$ , for $\theta>1$ I have to obtain: a)posterior distribution of $\theta$ given $Y$, where $Y = max(X1,...,X4)$. b)point estimate of $\theta$ based on median of posterior. **My work** a) Starting with pdf of $Y$, $f_{X(4)}(y)=\frac{4y^3}{\theta^4}$ ; using the formula for pdf of $X_{(n)}$ Calculating joint pdf $f(y,\theta)$ $f(y,\theta)=\frac{8y^3}{\theta^7}$ Posterior distribution $f(\theta|y) \propto \frac{1}{\theta^7}$ $f(\theta|y)=\frac{6}{\theta^7} \sim pareto(1,6)$ [By normalization] Could you please let me know whether my approach for posterior is correct. Also, for part (b), should I be obtaining the MLE for posterior distribution. I would appreciate your feedback, and pointers.


title: What is the difference between AIC() and extractAIC() in R?

text: The R documentation for either does not shed much light. All that I can get from this link is that using either one should be fine. What I do not get is why they are not equal. Fact: The stepwise regression function in R, `step()` uses `extractAIC()`. Interestingly, running a `lm()` model and a `glm()` 'null' model (only the intercept) on the 'mtcars' data set of R gives different results for `AIC` and `extractAIC()`.               > null.glm = glm(mtcars$mpg~1)     > null.lm = lm(mtcars$mpg~1)          > AIC(null.glm)     [1] 208.7555     > AIC(null.lm)     [1] 208.7555     > extractAIC(null.glm)     [1]   1.0000 208.7555     > extractAIC(null.lm)     [1]   1.0000 115.9434      It is weird, given that both the models above are the same, and `AIC()` gives the same results for both. Can anyone throw some light on the issue?


title: Random search for the optimal number of input features and optimal number of hidden layers for a MLP?

text: I've performed a random search in hypothesis space $$\{(c,h)| c \in U[1,256]; h\in U[1,100];c \in \mathrm{Z} \text{ and } h \in \mathrm{Z}\}$$ that defines the parameters of a standard multilayer perceptron (MLP) neural network. In each step of the random search, I draw two parameters $c$ and $h$. $c$ defines the number of input features and $h$ defines the number of hidden layer nodes. $c$ and $h$ are integers drawn from a uniform distribution defined above. I train a neural network defined by $(c,h)$ and calculate a misclassification rate and average squared error rate for each model. This is done with $10$-fold cross-validation to estimate the true error for each $(c,h)$. I therefore have an average misclassification rate and an average square error rate over the train sets and the left-out sets for each parameter pair. The question is, how do I chose the best pair of $(c,h)$ and is the method I use here sufficient? There is no reasonably clear point in the results as I'd have hoped. The results over the hypothesis space in the training data is: ![Training sample estimated errors](http://i.stack.imgur.com/ITpEf.png) The results over the hypothesis space in the hold-out data is, ![Hold out sample estimated errors](http://i.stack.imgur.com/kXlYE.png) This question relates to work I've done as part of my masters dissertation, and is related to the question here this


title: Clustering a binary matrix

text: I have a semi-small matrix of **binary features** of dimension 250k x 100. Each row is a user and the columns are binary "tags" of some user behavior e.g. "likes_cats".               user  1   2   3   4   5  ...     -------------------------     A     1   0   1   0   1     B     0   1   0   1   0     C     1   0   0   1   0      I would like to fit the users into 5-10 clusters and analyze the loadings to see if I can interpret groups of user behavior. There appears to be quite a few approaches to fitting clusters on binary data - what do we think might be the best strategy for this data?   * PCA   * Making a _Jaccard Similarity_ matrix, fitting a hierarchical cluster and then using the top "nodes".   * K-medians   * K-medoids   * Proximus?   * Agnes So far I've had some success with using hierarchical clustering but I'm really not sure it's the best way to go..               tags = read.csv("~/tags.csv")     d = dist(tags, method = "binary")     hc = hclust(d, method="ward")     plot(hc)     cluster.means = aggregate(tags,by=list(cutree(hc, k = 6)), mean)      ![enter image description here](http://i.stack.imgur.com/cyT6t.png)


title: How to perform t-test with huge samples?

text: I have two populations, One with N (no of observations)=38704 and other with N=1313662. These data sets have some 25 variables, all continuous. I took mean of each in each data set and computed the test statistic using the formula (t=mean difference/std error). The problem is of the degree of freedom. By formula of df=N1+N2-2 we'll have more freedom than the table can handle. Any suggestions on this? How to check the t statistic here. I know that the t-test is used for handling samples but what if we apply this on large samples.


title: Creating an index with objective and subjective variables

text: I have some objective and subjective variables. My aim is to create a statistical index which defines variables and interrelations of variables accurately. For instance, suppose that I am creating a spirituality index including demographic variables and also subjective variables like religious perceptions etc.. Considering the relations between all variables, how can I create an index value?What will be the mathematical presentation of this? Thanks in advance


title: Bayesian network fundementals

text: What is the formula of P(W) at this Bayesian network? ![enter image description here](http://i.stack.imgur.com/sZilQ.png) **PS 1:** Is that formula gives me P(W) or not: P(Cloudy) × P(Sprinkler|Cloudy) × P(Rain|Cloudy) × P(Wet Grass | Sprinkler, Rain) **PS 2:** How we can formula of P(R|S) ?


title: Are negative Detection Limit for concentrations consistent

text: We are monitoring Air Quality: essentially concentrations. We collect time series and we periodically calibrate analysers in order to control the process. Knowing the specifications or assessing it experimentally we have, for each analyser, some low threshold on the dynamic range (measurement extent) where data become meaningless because of the method variability, imperfections & limitations. Concentrations are physical real positive quantities. Therefore real concentration cannot be negative: whatever distribution it has, it should be positive. Analyser are real, and when measuring very low concentration, it can sometimes produce negative or null concentration due to variability of the process (we assume there is no bias in the method). In analytical chemistry question of Detection Limit often arises. It is a common practice to truncate part of the data series bellow some 'arbitrary' threshold (modelled or experimentally assessed). Of course this truncature bias data, which leads to higher expectation. But using meaningless data bias too. We know that dealing with DL are not the perfect solution, it is a filtering technique with advantages and drawbacks. My question is the following: When (1) null is an hard boundary that cannot be crossed (real concentration distribution is positive); (2) process variability close to the lower boundary of dynamic range of the instrument cannot be assumed: symmetric or concentration independent, is it reasonable to set a negative or null Detection Limit to truncate our data? If so, what hypothesis must hold?


title: Panel data: IV with different frequencies

text: I would like to estimate a standard logit (panel data): $\text{logit}(P(y_{i,t}=1))= \alpha + \beta_1 x^1_{i,t} + \beta_2 x^2_{i,t} +\epsilon_{i,t}$. The problem I am facing, however, is that $x^2$ is _yearly_ , whereas $y$ and $x^1$ are _monthly_. Aggregating $x^1$ to the yearly level is not an option. So I naively created 12 observations for each $x^2$ in order to make it monthly and be able to match it to every $y$ and $x^1$. Yet I realize that this might cause problems:   1. It may create serial correlation   2. It artificially inflates the number of observations and hence reduces the SE of $x^2$'s coefficient.   3. The SE of other coefs are not affected.   4. All coefficients remain unbiased. My **first question** is: is this _diagnosis_ of the potential problems correct? Are there other problems I did not think about? My **second question** is whether a solution to the problems outlined above is to use robust standard errors, clustering the observations by year. Thank you in advance,


title: How does one calculate confidence intervals on predictions generated by logit in Stata?

text: I'd like to make an assertion about whether individuals in my dataset exceed legal standards for commitment, which is largely determined by estimated risk of recidivism. I have estimated a logit model predicting recidivism. Here is sample code written trying to determine what percentage of individuals had predicted recidivism of at least 75% with 75% confidence:               .logit sexrecS [INDEPENDENT VARIABLES]     .predict sr1     .predict stdsr1, stdp     .gen byte sr1C1R1=0      .replace sr1C1R1=1 if 1/(1 + exp(ln((1-sr1)/sr1) - 0.67448975*stdsr1)) > .75     .sum sr1C1R1 if sr1!=. & sexrecS!=.      I later replace 0.67448975 with 1.281551566 to reflect a higher legal standard of proof, but the code generates more, not fewer, hits, so I know I've done something wrong.


title: Application of Hidden Markov Model to CRM

text: What are some of the applications of a HMM in the marketing field - specifically CRM and targeted marketing? For example, is it mainly for predicting an outcome given a sequence? Such as if web site visitor visits {Page 2, Page 54, page 23} what is the probability they will $<next>$ visit page 22? Or visit page 22 in that session? Is HMM a stand alone technique or normally coupled with other things?


title: Optimal importance sampling with ratio estimator

text: This is probably a stupid question, but here goes: so we want to approximate the following expectation: $$\mathbb{E}[h(x)] = \int h(x)\pi(x) dx$$ Where $h(x)$ is an arbitrary function and $\pi(x)$ is a distribution, also for simplicity, let's assume that we actually know the normalizing constant for $\pi(x)$. Of course we would like to sample from the optimal proposal distribution: $$g(x) = \frac{|h(x)|\pi(x)}{Z}$$ but of course this is not going to be a form we can sample from and we could not even compute the importance weights since we would need to know $Z$: $$ w(x) = \frac{Z \ \pi(x)}{|h(x)|\pi(x)}$$ But, if we assume $g(x)$ can be sampled from, can we use the ratio importance sampling estimator?: $$ \frac{\int w(x)h(x)g(x) dx}{\int w(x)g(x)dx}$$ To be clear, the estimator can also be written by letting $\{x^{(i)}\}_{i=1}^N$ be a set of samples from distribution with density $g(x)$ (magically). We let $$w(x^{(i)}) = \frac{1}{|h(x^{(i)})|}$$ making the final estimator: $$ \mathbb{E}[h(x)] \approx \frac{\sum_{i=1}^N w(x^{(i)})h(x^{(i)})}{\sum_{i=1}^N w(x^{(i)})} $$ So, is it true that the above estimator is asymptotically unbiased (consistent)? Or have I missed something? If it is indeed unbiased, could this then be used in conjunction with Monte Carlo approaches to sample from $g(x)$, since they can be used (in theory) to sample from any distribution known up to a normalizing constant. Edit: fixed a typo, also, I was able to prove this is consistant, so my new question is: is this a good idea? Are there any paper analyzing this? does it have a standard name?


title: Multilevel Meta-analysis with several treatments

text: I have data on a number of studies comparing response rates to various treatments of a disorder. There are six different treatments, and some studies only test a single treatment. For each treatment in each study, I have the response rate and the sample size (plus some other variables). Normally, if there were only two treatments, and each study compared the two, I might use something like $$D_i=\beta_0+\beta_1\cdot SS_i+u_i+\epsilon_i$$ where $D_i$ is the difference between the logits of the response rates in study $i$, $\epsilon_i$ is the corresponding sampling error, $SS_i$ is the sample size of study $i$, and $u_i$ is the residual error term. Is it possible to sensibly analyze my data in a similar way? Are there keywords / references I should be aware of? For example, I'd like to use something like $${\rm logit}(RR_{ij})=\beta_0+\beta_1\cdot SS_{ij}+\beta_2\cdot T_j+u_i+\epsilon_{ij}$$ where $i$ is the study and $j$ is the treatment (1 to 6), $T_j$ is an indicator variable for the treatment (alternatively, could be turned into a second random effect), and $u_i$ and $\epsilon_{ij}$ are as above.


title: Measuring 'synchrony' with time series correlations

text: Question from a stats novice and StackExchange newb, an anthropologist with a computer science background. I'm looking for an appropriate statistical measure for 'synchrony' as in the following situation: I have a panoramic video of several people sitting in chairs around the camera having a conversation for an hour. I split the video up by person and generated time series data for each individual, measuring how much they move, gesture, etc. over time (i.e. number of pixels different between successive video frames). So I have a graph of how still or animated a person is at each moment over an hour (using moving average to smooth things out). I want to measure the 'synchrony' of the group at each moment, i.e. generate another time series that shows how 'synched up' their movements are at any given moment. What kind of statistical measure would be appropriate for this? Simpler/intuitive the better. Some assumptions:   * If only one person is moving, synchrony should be near zero.   * If everyone is moving, synchrony should be near max. Note that this will be used primarily for _visual_ inspection and analysis, to highlight those parts of the video where people are synchronous. -Daniel


title: How to interpret output of Match() function in R (for propensity score matching)

text: I'm trying to use the Match() function from the Matching package in R to do a propensity score analysis. My outcome of interest is a binary variable (0/1). My treatment is also a binary variable (0/1). In addition, I have a number of other variables that I want to control for in this analysis. First, I fit a logistic regression to define a propensity score for the treatment:               glm1 = glm(Treatment ~ variable1 + variable2 + variable3 + ...,                 data=dataset, family="bionomial")      Then, I used the Match function to estimate the average treatment effect on the treated:               rr1 = Match(Y = Outcome, Tr = Treatment, X = glm1$fitted)      Finally, I called for a summary:               summary(rr1)      My question is how to interpret the output. I get:               Estimate... -0.349,     AI SE... 0.124,     T-stat... -2.827,     p.val... 0.005      What does this mean? In particular, what is Estimate? The documentation says it's "The estimated average causal effect." But what are the units? Can I interpret this to mean that the treatment reduced the outcome by a relative 35%? Or by an absolute 0.35? Or do I need to exponentiate? Any help on the interpretation would be much appreciated!


title: testing logistic regression coefficients using $t$ and residual deviance degrees of freedom

text: **Summary:** Is there any statistical theory to support the use of the $t$-distribution (with degrees of freedom based on the residual deviance) for tests of logistic regression coefficients, rather than the standard normal distribution? * * * Some time ago I discovered that when fitting a logistic regression model in SAS PROC GLIMMIX, under the default settings, the logistic regression coefficients are tested using a $t$ distribution rather than the standard normal distribution.$^1$ That is, GLIMMIX reports a column with the ratio $\hat{\beta}_1/\sqrt{\text{var}(\hat{\beta}_1)}$ (which I will call $z$ in the rest of this question), but also reports a "degrees of freedom" column, as well as a $p$-value based on assuming a $t$ distribution for $z$ with degrees of freedom based on the residual deviance -- that is, degrees of freedom = total number of observations minus number of parameters. At the bottom of this question I provide some code and output in R and SAS for demonstration and comparison.$^2$ This confused me, since I thought that for generalized linear models such as logistic regression, there was no statistical theory to support the use of the $t$-distribution in this case. Instead I thought what we knew about this case was that   * $z$ is "approximately" normally distributed;   * this approximation might be poor for small sample sizes;   * nevertheless it _cannot_ be assumed that $z$ has a $t$ distribution like we can assume in the case of normal regression. Now, on an intuitive level, it does seem reasonable to me that if $z$ is approximately normally distributed, it might in fact have some distribution that is basically "$t$-like", even if it is not exactly $t$. So the use of the $t$ distribution here does not seem crazy. But what I want to know is the following:   1. Is there in fact statistical theory showing that $z$ really does follow a $t$ distribution in the case of logistic regression and/or other generalized linear models?   2. If there is no such theory, are there at least papers out there showing that assuming a $t$ distribution in this way works as well as, or maybe even better than, assuming a normal distribution? More generally, is there any actual support for what GLIMMIX is doing here other than the intuition that it is probably basically sensible? R code:               summary(glm(y ~ x, data=dat, family=binomial))      R output:               Call:     glm(formula = y ~ x, family = binomial, data = dat)          Deviance Residuals:         Min      1Q  Median      3Q     Max       -1.352  -1.243   1.025   1.068   1.156            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  0.22800    0.06725   3.390 0.000698 ***     x           -0.17966    0.10841  -1.657 0.097462 .       ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1235.6  on 899  degrees of freedom     Residual deviance: 1232.9  on 898  degrees of freedom     AIC: 1236.9          Number of Fisher Scoring iterations: 4      SAS code:               proc glimmix data=logitDat;         model y(event='1') = x / dist=binomial solution;     run;      SAS output (edited/abbreviated):               The GLIMMIX Procedure                         Fit Statistics          -2 Log Likelihood            1232.87     AIC  (smaller is better)     1236.87     AICC (smaller is better)     1236.88     BIC  (smaller is better)     1246.47     CAIC (smaller is better)     1248.47     HQIC (smaller is better)     1240.54     Pearson Chi-Square            900.08     Pearson Chi-Square / DF         1.00                                      Parameter Estimates                                   Standard     Effect       Estimate       Error       DF    t Value    Pr > |t|          Intercept      0.2280     0.06725      898       3.39      0.0007     x             -0.1797      0.1084      898      -1.66      0.0978      $^1$Actually I first noticed this about _mixed-effects_ logistic regression models in PROC GLIMMIX, and later discovered that GLIMMIX also does this with "vanilla" logistic regression. $^2$I do understand that in the example shown below, with 900 observations, the distinction here probably makes no practical difference. That is not really my point. This is just data that I quickly made up and chose 900 because it is a handsome number. However I do wonder a little about the practical differences with small sample sizes, e.g. $n$ < 30.


title: Error bars calculation in a split-plot design

text: I have the following experimental design :   Sampling is carried out on three random sites (Factor Site is coded : Site1, Site2, Site3).   Within each site, sampling is carried out in a dry plot and in a wet plot (Factor Hydrology is coded : Dry/Wet)   Within each of the previous plot, sampling is carried out at two Depths (D1, D2) in triplicate.   There are a total of n = 36 samples = 3 Sites * 2 Hydrology * 2 Depths * 3 Replicates. I use the lme function in R (lnme package) as follows :               Site<-as.factor(rep(c("Site1","Site2","Site3"),each=12))     Hydrology<-as.factor(rep(rep(c("Dry","Wet"),each=6),3))     Depth<-as.factor(rep(rep(c("D1","D2"),each=3),6))     Variable<-rnorm(36)               mydata<-data.frame(Site,Hydrology,Depth,Variable)          mod1<-lme(Variable~ Hydrology*Depth, data=mydata, random=~1|Site/Hydrology/Depth)      I want to visualize the impact of Hydrology and Depth on the dependent variable without displaying the data for each Site :               interaction.HD<-interaction(mydata$Hydrology,mydata$Depth)     calculated_mean<-tapply(mydata$Variable,interaction.HD,mean)     mean_V=as.table(matrix(calculated_mean,nrow=2))     rownames(mean_V)=c("Wet","Dry")     colnames(mean_V)=c("D1","D2")          library(gplots)     graph=barplot2(mean_V, beside = TRUE)     legend("center",c("Dry","Wet"),pch=22,pt.bg = c("red","yellow"))      Each bar on the graph represents the mean of 9 values (3 replicates per Site for each combination of Hydrology and Depth). I would like to plot error-bars on each mean. My question is the following : What is the most correct way to calculate it? I think that it is incorrect to calculate the Standard-Deviation or the Standard-Error on the 9 values used to calculate each mean. I think it would be more correct to calculate SD or SE on 3 values (each value corresponding to the mean per site for each combination of Hydrology and Depth). Any advice on this question would be much appreciated. Do you think it is correct to use the method described in Cousineau (2005) where the between-subject variability is removed before the calculation of error bars? link to the reference : http://www.tqmp.org/Content/vol01-1/p042/p042.pdf


title: What kind of distribution is it?

text: I was trying to help my colleague with fitting a distribution curve to some empiric data. However, I haven't succeeded since it's my first experience in approximating an empiric distribution with a theoretical one. The mean value of $Y$ is 0.32688. $\sigma^2_Y$ is 0.076191. Here's a histogram: ![Y histogram](http://i.stack.imgur.com/LDB4m.png) And an empirical cdf: ![Y empirical cdf](http://i.stack.imgur.com/OxS1q.png) There are 1024 samples in the data. Obviously the distribution is asymmetric. It doesn't look like $\chi^2$ since $\hat{\mu}_Y >> 2\hat{\sigma}_Y^2$ I tried approximating it with beta distribution (picked various $\alpha, \beta$ parameters) together with normalizing Y by some real number in $[0.5; 3]$. It didn't help, Pearson chi-squared statistic was always, you know, like 300, 700 or even more. I tried squaring Y and the histogram of $Y^2$ looks like this: ![enter image description here](http://i.stack.imgur.com/QyvOK.png) I tried approximating this with exponential distribution (and normalizing Y with a real number) - but P-value was always like zero all the time. The data itself can be found here: https://docs.google.com/spreadsheet/ccc?key=0AqeuZRXRoMfRdF9UQzFEZVdOTmh0ZzNSaFRhRndHWmc Does anyone have any ideas on what distribution this data could have?


title: A name for this distributional condition?

text: I have come across a needed condition on a continuous probability distribution defined over $[0, \infty]$ and wonder whether it has a name. For a distribution with CDF $F$ and pdf $f$ I need that the quantity: $$\phi(x) \equiv \frac{f(x)}{F(x)+xf(x)}$$ be monotonically non-increasing. Putting the condition in another form (by taking derivatives), the requirement is that for all $x \in [0,\infty]$ such that $f'(x) > 0$: $$f(x) \geq \sqrt{\frac{F(x) f'(x)}{2}}$$ This seems to be a commonly satisfied property, so does it have a name? It is related to, but distinct from a monotone hazard rate condition.


title: "Standard Deviation" of non-negative data

text: I've been asked to analyse the distribution of a set of data, essentially a single column of random samples of a physical parameter which cannot be negative. The standard deviation that I usually calculate assumes that data are normally distributed, which in this case cannot be true. Since I'm used to being asked to provide the standard deviation and mean of a dataset as a measure of the errors, I'd like to know whether there is an easily calculated parameter for a more appropriate distribution that would fulfil the same purpose? I apologise if this question is ill-posed; I am not a statistician by training. With many thanks, Loruschorus


title: Probabilities in a Markov Model

text: I am reading a paper on Markov Models and I am trying to figure out how to compute the probabilities for the $\alpha$-pass. I am given an $N\times N$ matrix $A$, that has the probabilities of transitioning from one state to to another. An $N \times M$ matrix $B$, that has the probabilities of a state, given an observation $M$. In my model I have $N=2$ and $M=3$ I have an initial state sequence $\pi = (\pi_{x_0},\pi_{x_1})$ I also have a set of observations $O=(O_0,O_1,O_2,O_3)$. Let $X=(x_0,x_1,x_2,x_3)$ be a state sequence. I know that $P(O,X|\lambda)=b_{x_0}(O_0)b_{x_1}(O_1)b_{x_2}(O_2)b_{x_3}(O_3)$ $P(X|\lambda)=\pi_{x_0}a_{x_0,x_1}a_{x_1,x_2}a_{x_2,x_3}$ I also know that: $P(O,X|\lambda) = \pi_{x_0}b_{x_0}(O_0)a_{x_0,x_1}b_{x_1}(O_1)a_{x_1,x_2}b_{x_2}(O_2)a_{x_2,x_3}b_{x_3}(O_3)$ I am interested in the probability: $\alpha_t(i) = P(O_0,O_1,O_2,O_3,x_t=q_i|\lambda)$, where $t=0,1,2,3$ and $i=0,1$, $q_i$ is a state. I have that: $\alpha_0(i) = \pi_i b_i(O_0)$ for $i=0,1$ My question is how do I calculate $\alpha_1(i), \alpha_2(i), \alpha_3(i)$?


title: Document classification problem

text: > Assume we have $L$ labelled documents, and $U$ unlabeled ones, where all the > documents from class $k$ were generated from a multinomial or Naive Bayes > distribution with parameter $\theta_k$, and there are a total of $K$ > classes. Let each point $x_n$ be a $D$ dimensional vector. > > a) What is the maximum likelihood estimator for $\theta_k$ using only the > labeled data? > > b) What is the mean squared error of this estimator? > > c) We can construct an estimator for $\theta_1, \ldots, \theta_K$ that uses > the unlabeled as well as labeled data by maximizing the likelihood of the > observed data. Write down the log-likelihood function in this case. > > d) Under what conditions would the MLE in case (a) converge to the true > parameter? Can these conditions be relaxed if we also use unlabeled data as > in (c)? I am not sure as to how to solve each of these. For part (a), is this the likelihood $$ l(\theta_{k,d}) = \prod_{n=1}^N p(y_n)^{I(y_n = k)} \left[ \prod_{d=1}^D p(x_n^d | y_n = k) \right]^{I(y_n = k)} $$ that I should maximize with respect to $\theta_{k,d}$? If so, is the estimator $$ \theta_{k,d} = \frac{N_k}{N} $$ where $N_k$ is the number of data points labeled $y_n = k$? Can anyone verify this? And can you give me suggestions as to how to solve the remaining parts? Thanks


title: Survival analysis w/ repeated observations of experimental and control units within periods of differing size & trends over the periods

text: I have an experiment called Fear Conditioning Acquisition Protocol. We expose mice to subsequent white sounds followed by short electric shocks through the floor. Our aim is to let them learn to associate the sound to the shock. We assess the degree of learning by measuring the freezing of the mouse during the sound exposition. To do the comparison between the subsequent phases of the experiments we use the average of the freezing of each mouse in every phase. This is the result graph: ![Exposed group](http://i.stack.imgur.com/Zm9rr.png) Then we have another set of mice in which no shock follows the sound (control group) here's the results: ![Control group](http://i.stack.imgur.com/MxHVq.png) Values are percentage of freezing on phase duration. My questions are:   1. which test should I use to compare the average freezing between each phase to assess if it's significantly different? I used ANOVA with "Bonferroni's Multiple Comparison Test".    2. graphically the level of freezing is clearly increasing, is it important then to question about the significance of the difference among the freezing levels?    3. how to formally compare the freezing in the exposed and the control group to show that is significantly different?   4. Generate a value that represent the freezing of a mouse to correlate it with the expression of certain proteins. What to use? total freezing over total sound exposition? average freezing?


title: How to report correlation between two variables across multiple groups in a table?

text: Following on from my post here, I have three statistically significant correlations between Variables A and B for three age groups. They are:   * Group 1 (Less than 18 years) r1 = .74 p < .001 n = 99   * Group 2 (18 to 65 years) r2 = .78 p < .001 n = 96   * Group 3 (More than 65 years) r3 = .75 p < .001 n = 97 Note: I have done Šidák correction on the alpha threshold of 0.05. However, I am most interested in comparing the correlation coefficients and I have used the calculator here to do the analysis which gives a $p$-value of .794. This is not statistically significant at the new alpha threshold of .0006 (so $H_0$ is supported i.e. no difference between the groups, though there is a significant relationship between A and B for each group). **Question 1: How do I report these significant and this non significant findings in a single table?** I have columns labelled (in this order):   * Group   * r   * r squared   * p value   * sample size **Question 2: Should I also include the confidence interval. How does this add to the information already provided, especially when I am noting the p value?** One of the contributors below suggests that there may be an alternative way of addressing the problem above (which I can only presume to be a better/robust method than Pearson's r that I have used. **Question 3: Is there an alternative analysis I can perform for the problem above (noting that I am focusing on relationships and the equality of the relationships across demographic characteristics)?** (I think by now I really understand correlation is not causation; I am not interested in causation!) Some insights into why I am asking these questions, especially Q3, are in response to another question here.


title: Robust regression: possible with ordered data?

text: I have two regressions to perform - one with a metric DV (-3 to 3), the other with an ordered DV (0,1,2,3). Neither normal distribution nor homoscedasticity is given. I have a two questions:   1. Some sources say robust regression take care of both lack of normal distribution and heteroscedasticity, while others say only of normal distribution. What is true?    2. Are there ways of using robust regressions with ordered data, or is that only possible for metric DVs?  I have SPSS, R, and MPLUS available.


title: Signal extraction problem: conditional expectation of one item in sum of independent normal RVs

text: I have what I'm afraid is a simple stats problem that is stumping me. I have two random variables, X and Y, independently normally distributed:                X ~ N(0, sigmaX)      Y ~ N(0, sigmaY).      I observe the sum of these two variables, Z = X+Y, and want to develop a conditional expectation on X given the sum. A colleague said, "ah, yes, classic signal-extraction problem. Solution is:"                E[X|X+Y] = (X + Y) * sigmaX / (sigmaX + sigmaY)      This looked about right so I thanked him and figured I work it out at home. It appears I'm a little rusty here, though. I can give verbal reasoning why this would be true but can't write down the math. What is the mathematical reason this is true? Thanks all!


title: Correlating time stamps

text: Let's say I am analyzing behavioral patterns over the course of an hour. I have recorded three different behaviors and the time stamps (start end) they occurred at. Something like:               yawning       stretching    whispering     2:21-2:22     3:31-3:33     1:21-1:30     3:42-3:45     8:23-8:59     9:27-9:33     9:20-925      9:34-9:44     14:04-14:07     14:45-14:32   15:01-15:06   18:00-18:22     .     .     .     45:40-45-43   45:23-45:30   44:19-44:44      Is there a statistical method for determining if certain behaviors correlate or cluster around certain time periods/to each other? For instance maybe I want to know if these three (or just 2) behaviors are found in close proximity to one another or maybe I want to know if these which behaviors are not in close proximity to each other. Which of the three behaviors tend to cluster together? I don't even know what field of stats I'm looking at with this.


title: Mismatch between significant variables from logistic regression and tree splits in R

text: I'm studying a data set in R using both regression trees (tree and rpart functions) and logistic regression. I'm finding explanatory variables in the regression which are significant, but when I fit a tree those variables are not used as splits. What does this imply about the predictive ability of the results?


title: What is the best way to visualize a single numeric variable as a heatmap?

text: I'm a hobbyist programmer whose friend recently took a business trip overseas. He's polled our mutual friends for bets on the size of his email inbox when he returns. I'd like to visualize this as a one-dimensional heatmap where color reflects a kind of density of guesses for that value. I understand I could create a histogram and find a way to translate its values into color values. However, the histogram bins would create janky transitions in color. Is there a way to create an appropriate, continuous function for prettier output?


title: How to run and plot simulation in R of sum of 20 random variables?

text: I would like to run a simulation where the sum of 20 random variables are generated 1000 times and plotted in a histogram. I generate the random numbers using runif(20,0,1) then sum these random numbers using the sum function, but I was wondering if there is an easier method to do this 1000 then to write a for loop.


title: Does there exist a model fit statistic (like AIC or BIC) that can be used for absolute instead of just relative comparisons?

text: I'm not that familiar with this literature, so please forgive me if this is an obvious question. Since AIC and BIC depend on maximizing the likelihood, it seems that they can only be used to make relative comparisons between a set of models attempting to fit a given data-set. According to my understanding, it wouldn't make sense to calculate the AIC for Model A on data-set 1, calculate the AIC for Model B on data-set 2, and then compare the two AIC values and judge that (for example) Model A fits data-set 1 better than Model B fits data-set 2. Or perhaps I am mistaken and that is a reasonable thing to do. Please let me know. My question is this: does there exist a model fit statistic that can be used for absolute instead of just relative comparisons? For linear models, something like $R^2$ would work; it has a defined range and discipline specific ideas about what is a "good" value. I'm looking for something more general and thought that I could start by pinging the experts here. I'm sure that someone has thought of this kind of thing before, but I don't quite know the right terms to make a productive search on Google Scholar. Any help would be appreciated.