title: Issue in training Hopfield network and convergence problem

text: I am learning how to use Hopfield Neural network as a context addressable memory. The objective is to obtain a fixed point of the network which indicates an equilibrium state. This state vector remains unchanged for successive iterations and is called a fixed point. When the network reaches the fixed point, we say that a decision or goal has been reached. As can be seen from below illustration, by recursive application of the formula, X2 is the fixed point.               WeightMatrix           = [0.0 0.0 -1 0 1;                               0.0 0.0 0 -1 0;                               0.0 -1 0.0 0 -1;                              -1 1.0 0.0 0.0 0.0;                               0.0 0.0 0.0 1 0.0];          X1 = (1 0 0 0 0)   %Initial training example     X1*W = [0,0,-1,0,1];          X2 = f(X1*W) = [1,0,0,0,1]     X2*W = [0,0,-1,1,1];  X3 = f(X2*W) = [1,0,0,1,1]     X3*W = [-1,1,-1,1,1]; X4 = f(X3*W) = [1,1,0,1,1]     X4*W = [-1,1,-1,0,1]; X5 = f(X4*W) = [1,1,0,0,1]     X5*W = [0,0,-1,0,1];  X6 = f(X5*W) = [1,0,0,0,1] = X6 = **X2**      Problem : The above example is based on the paper Download link titled "Application study in decision support with fuzzy cognitive map" explains how to train a [Fuzzy Cognitive Map][2] with this example in Section 2.2. When I simulated the example, using the same logistic sigmoid thresholding function               logistic sigmoid function f(X1*W) = 1/(1+exp(-X1*W).       , I am not getting the same output after passing the result into the thresholding function as mentioned in the paper. As a result, the program/network is iterating millions of time and not converging to a fixed point. By passing the result of X1*W into f, I am getting real valued numbers for X2 and not binary !! Same problem for the rest of the iterations. Is there something wrong in my understanding of how convergence/memory recall is performed or is my code incorrect? Please help CODE               Training1 = [1,0,0,0,0];          lambda =1;      t = 1;     X(t,:) = Training1;      err = 1;       while(err~=0)       Out = X(t,:)*WeightMatrix;       temp  =  1./(1+exp(-lambda.*Out));         err = ((Out - temp)*(Out - temp).')/numel(temp);     t = t+1      X(t,:) = temp;      end


title: Questions on linear regression

text: I had a few true / false questions on a practice test that I would like to discuss if possible.   * A value $R^2$ close to 1 indicates the linear regression is a good fit to data Yes, but I am not sure how to reason it other than intuitively. More precisely, why do we decide to define $R^2$ as 1-(sum of residuals/total sum of squares)?   * The estimate of the error variance $s^2$ is a random variable I believe it is because how the error varies is normally distributed if I remember correctly. If anyone could elaborate it would be much appreciated.   * $∑(Y_{1,i}−Y_{2,i})^2 = 0$, where $Y_1$ is the predicted value and $Y_2$ is the actual value. Since the sum of residuals is zero and since this is equivalent to the sum of residuals it must be zero.


title: corpus linguistics statistics

text: ![enter image description here](http://i.stack.imgur.com/Dq3fz.jpg) I am in doubt as to how calculate observed relative frequences in corpus linguistics. That's how I did it: I multiplied the number of words by thousand and then I divided it all by the size of the corpus, for example, 81*1000/130621= 0.62. I did the same thing with the other corpus of a different size: 49*1000/33756=1.4. Does it seem to you that corpora is normalized/standartized? As for me, it does not seem to be because of the higher ratio in the second corpus (popular scientifc articles'). Thanks


title: What options are there to combine different distance functions?

text: I am currently working with feature-vectors that are made up of continuous attributes, so I can use the euclidean distance for things like KNN- classification and clustering. Now I want to add a nominal attribute that has a special distance function defined. What options do I have of combining these distance functions, so I still obtain one distance for two vectors?


title: Confidence interval for ratio in timeseries

text: I have a time series X(t). Each X(t) has three possible outcomes : A,B,C. I am interested in the ratio of A,B,C to the total. Assuming N is the total number of data points I have gathered for X(t),   _How can I compute the confidence levels for A/N, B/N and C/N when the X(t) are "intuitively" not independent ?_ For exemple : X(t) is an indication of whether a car is : moving (speed>0), stopped (speed=0) or its motor is off. The data I gathered is a time slice pertaining to a car. To me, those categories are not independent because when the car is moving at X(t) is is very likely to be still moving at X(t). [Am I correct?]


title: How to calculate minF' with corrected degrees of freedom?

text: If you're working with crossed random effects (a sample of words and a sample of people), you should report a minF' rather than just your F1 and F2 values according to Clark (1973), Raaijmakers et al. (1999), ... I've found this site which calculates the minF' for me. However, this only works if my degrees of freedom for F1 and F2 are the same. Since the assumption of sphericity was violated in my F1 analysis, I used a Greenhouse- Geisser correction, which resulted in different degrees of freedom. How can I calculate a minF' in this case? Or do the degrees of freedom not play a role, given an F-statistic? This is possible, since the original formula only worked with the F-statistics. Is there a site or a piece of software (similar to the one I linked) that can do this?


title: Is the power for a Kappa test the same as underlying z-test?

text: The Kappa ($\kappa$) test is a Z-test kind of test. If I am not very wrong, to compute the $\kappa$ test, we can just estimate the appropriate variance $\hat {var}(\hat\kappa)$ for the kappa statistic $\hat\kappa$ and then feed it to a z-test by taking $\mu$ = $\hat\kappa$ and $\sigma^2$ = $var(\hat\kappa)$. To compute the power of a z-test one can use the relation $1 - \beta = \phi(Z_{a} - \sqrt n * (\mu-\mu_0)/\sigma)$ Would the power of this underlying z-test be also the power of the original $\kappa$ test? If not, why?


title: Do third order asymptotics exist?

text: Most asymptotic results in statistics prove that as $n \rightarrow \infty$ an estimator (such as the MLE) converges to a normal distribution based on a second-order taylor expansion of the likelihood function. I believe there's a similar result in Bayesian literature, the "Bayesian Central Limit Theorem", which shows that the posterior converges asymptotically to a normal as $n \rightarrow \infty$ My question is - does the distribution converge to something "before" it becomes normal, based on the third term in the Taylor series? Or is this not possible to do in general?


title: Where do the full conditionals come from in Gibbs sampling?

text: MCMC algorithms like Metropolis-Hastings and Gibbs sampling are ways of sampling from the joint posterior distributions. I think I understand and can implement metropolis-hasting pretty easily--you simply choose starting points somehow, and 'walk the parameter space' randomly, guided by the posterior density and proposal density. Gibbs sampling seems very similar but more efficient since it updates only one parameter at a time, while holding the others constant, effectively walking the space in an orthogonal fashion. In order to do this, you need the full conditional of each parameter in analytical from*. But where do these full conditionals come from? $$ P(x_1 | x_2,\ \ldots,\ x_n) = \frac{P(x_1,\ \ldots,\ x_n)}{P(x_2,\ \ldots,\ x_n)} $$ To get the denominator you need to marginalize the joint over $x_1$. That seems like a whole lot of work to do analytically if there are many parameters, and might not be tractable if the joint distribution isn't very 'nice'. I realize that if you use conjugacy throughout the model, the full conditionals may be easy, but there's got to be a better way for more general situations. All the examples of Gibbs sampling I've seen online use toy examples (like sampling from a multivariate normal, where the conditionals are just normals themselves), and seem to dodge this issue. * Or do you need the full conditionals in analytical form at all? How do programs like winBUGS do it?


title: Can conditional distributions be considered as a parameterized family of distributions?

text: Given two real-valued random variables $X$ and $Y$, are the conditional distributions of $Y$ given $X$ taking different values, $ \\{ p(y|x), \forall x\\}$, considered as a parameterized family of distributions, i.e. a parameterized model (in statistics)?


title: Hypothesis testing simple null vs. complex alternative

text: This is a simple hypotheses testing question between a simple null and a complex alternative. I want to test the null hypothesis $H_0 : \mu = \mu_0$ against the alternative $H_1 : \mu = \mu_0 \pm \epsilon =\mu_1$ where $\epsilon$ is a fixed number where $\mu_0+\epsilon$ and $\mu_0-\epsilon$ are both equally likely. So can I say that the most powerful test is the one where we test for $\mu_0>\mu_1$ and $\mu_0<\mu_1$?


title: Using day length and rate of change of day length in linear regression

text: I have a time series in which samples were collected approximately monthly over ten years. I have wanted to identify whether a number of bacterial species were seasonally variable and whether they increased or decreased over time. I had been modelling species abundance with linear regression according to the equation S = a + b*ED +c*DL + d*DDL. Where S is the abundance of a species. ED: Number of days elapsed since the beginning of the study. DL: Day length on the day the samples were collected DDL: The rate of change in day length, measured here as the day length minus the day length measured thirty days previous. I then looked at the P values associated with each of these coefficients and the coefficients themselves. If ED was statistically significant, I said there was a long term triend, and if ED's coefficient was positive, that species abundance was increasing. If DL or DDL's coefficients were statistically significant, I said the species' abundance was seasonal. If DL had a statistically significant coefficient, I said the species was most abundant in the summer; if negative, most abundant in the winter. Same pattern for DDL, positive, most abundant in spring, negative most in fall. A reviewer tells me that "The regression analysis may not be appropriate because day length and day length change are autocorrelated and also covarying." Is this criticism valid? If valid, is there anything I can salvage from this kind of analysis? Is there a better approach? If not, why isn't it valid? I feel like this approach is analogous to building trigonometric models to describe a data set eg y = a + b*t + c*sin(2pi*t/L) + d*cos(2pi*t/L), which people apply to model some time series. Can I just replace the sign and cosine terms with DL and DDL (which are sinusoidal and DDL is basically the derivative of DL) or am I breaking some rule there? Thanks!


title: Q-Q plot and KS test

text: I have a question regarding the QQ plot and KS test. I have a sample which is some positions on genome, and I calculated the distance between these positions. And I generated the null distribution by randomly draw positions from genome and calculated the distance between them in the same way. (1000 times) I want to test if the distribution of the distance of my sample are the same with the distribution of distance of randomly selected positions. First I used the QQ plot to test it, and I get the QQ plot looks like this: ![QQ plot on log10 scale](http://i.stack.imgur.com/aiI2h.jpg) The QQ plot is skewed, then I want to use KS test to see if they are from the same distribution. However, I could not get a significant P-value from using ks.test() function in R. The P-value is around 0.09... To be clear, this is the histogram of null distribution I generated by simulation: ![enter image description here](http://i.stack.imgur.com/CTUAV.jpg) Does anyone have idea why I couldn't get significant result from KS-test? Or if there is any alternative method I could use to test the difference between this two distributions? Thanks so much!


title: Details on kriging and variograms

text: I was wondering what are the steps to perform factoring kriging. I know it is used to extract specific components that may have a physical meaning and this is achieved by kriging. However, I have some questions regarding the specifics of such procedure. By the way, the following description of the problem belongs to Wackernagel's excellent book Multivariate Geostatistics. Suppose we have a random variable $Z({\bf{x}})$ that can be decomposed into some stationary components $Z_{u}({\bf{x}_{\alpha}})$ and an intrinsic component $Z_{S}({\bf{x_{\alpha}}})$. We are interested in kriging this intrinsic component: $Z_{S}^{*}({\bf{x_{0}}}) = \sum_{\alpha=1}^{n}w_{\alpha}^{S}Z({\bf{x_{\alpha}}})$ Therefore, estimating the variance $\sigma_{E}^{2}$ we get: \begin{array}{} \sigma_{E}^{2} &= var(Z_{S}^{*}({\bf{x_{0}}})-Z_{S}({\bf{x_{0}}}))\\\ & = \sum_{u=0}^{S-1}C^{u}({\bf{x_0}}-{\bf{x_0}}) - \gamma^{S}({\bf{x_0}}-{\bf{x_0}}) - \sum_{\alpha=1}^{n}\sum_{\beta=1}^{n}w_{\alpha}^{S}w_{\beta}^{S}\gamma({\bf{x_\alpha}}-{\bf{x_\beta}})\\\ &+2\sum_{\alpha=1}^{n}w_{\alpha}^{S}\gamma^{S}({\bf{x_\alpha}}-{\bf{x_0}}) \end{array} and minimizing (with respect to the weights) we obtain: $$\sum_{\beta=1}^{n} w_{\beta}^{S}\gamma({\bf{x_\alpha}}-{\bf{x_\beta}}) + \mu_{S} = \gamma^{S}({\bf{x_\alpha}}-{\bf{x_0}}) \thinspace \thinspace \text{for $\alpha$=1,...n}$$ $$\sum_{\beta=1}^{n}w_{\beta}^{S}=1$$ where $\mu_{S}$ is the Lagrange multiplier due to the constraint on the weights. As I understand it, the idea usually is to get the weights $w_{\beta}$ that minimize the variance $\sigma_{E}^{2}$ in order to estimate new values ${\bf{x_{0}}}$. In this case, the values to be estimated are the intrinsic components $Z_{S}^{*}({\bf{x_{0}}})$. However, I'm not sure how to fit from data the variograms needed for these equations. Furthermore, in factoring kriging the variogram or covariance function is decomposed into several variograms that correspond to the components of $Z({\bf{x}})$. Assuming I know the variograms, should I use something like $\gamma_{1}({\bf{x_\alpha}}-{\bf{x_\beta}}) + \gamma_{2}({\bf{x_\alpha}}-{\bf{x_\beta}}) ... $ in place of $\gamma({\bf{x_\alpha}}-{\bf{x_\beta}})$ in the kriging equations?


title: Correlation between predictor variables and time in Survival Analysis?

text: I am using survival analysis to model time to an event. I would like to explore the effect of a continuous predictor variable on the hazard rate. The continuous predictor ranges in value from 0-8.5. This variable is fixed (1 value) for each observation, i.e. not time dependent. The data is also characterized by left-truncation, i.e. some patients enter into the study after time zero. The issue I am having is that the independent variable is not so independent. The upper portion of its range is only ever associated with small values of time. I have convinced myself that this dependence issue precludes me from using this type of analysis. Can anyone point to reference that may convince me otherwise.


title: What can be inferred from a 95% confidence interval on a correlation coefficient?

text: If the only information you have for a pearson's correlation is the 95% confidence interval, what can you infer from that data? For example, if you had a correlation coefficient of (0.24;0.78) what would be the best inference to make? I don't have a strong background in stats so if someone could explain it without lots of equations that would be preferable, thanks!


title: Visualizing many left-skewed distributions

text: I have a series of left-skewed/heavy tailed distributions that I would like to show. There are 42 distributions across three factors (labeled as `A`, `B` and `C` below). Also, the variation is shrinking across factor `B`. The issue I have is that the distributions are hard to differentiate across the scale of the outcome (a ratio or fold-change): ![enter image description here](http://i.stack.imgur.com/qImU4.png) Logging the data seems to over-emphasizes the left skewness and moves more samples into the tails (creating a mash of outlier points): ![enter image description here](http://i.stack.imgur.com/RBVJM.png) Does anyone have suggestions on other techniques for visualizing these data? Thanks.


title: Is a nested log-binomial GLM appropriate for this?

text: Dependent variable: Death (0 or 1)   Independent variables: Dose (0, 10 or 20), Timing(Early or Late), PreviousDVT(0 or 1) Ignoring the Timing variable for the moment, the model I had originally considered doing was a log-binomial GLM with the formula Death ~ Dose * PreviousDVT There (should) be all combinations of Dose X PreviousDVT available in the data set (Dose 0 with/without previous DVT, dose 10 with/without previous DVT and dose 20 with/without previous DVT). Here's where I'm running in to difficulty. I also want to know whether giving this drug within 24 hours (early) affects the outcome as compared with giving it after 24 hours (late). If the dose is 0, then Timing simply does not apply. I presume I need some sort of nested design (Death ~ PreviousDVT * Dose/Timing), but it's been a while and I'm also not sure how to ensure that when I look at Timing I don't get spurious results due to all the Dose=0 cases contributing. To summarise:   (1) I have a variable where one level (Dose=0) means that the levels in another vairable (Timing) are meaningless. Do I need another level in Timing, i.e. "not applicable"?   (2) How do I model this problem in R? It's been a while since I've done nested desgins (presuming this is one) - is it Death ~ PreviousDVT * Dose/Timing? Thanks very much!


title: Applying autoencoders for dimensionality reduction in audio: Why does this create a low-pass effect?

text: I've been playing around with framing audio data and training a single-layer autoencoder to find a dimensionality-reduced form (say 128-sample frames to 32-dimension frames). When I test the audio after going through the encode/decode layers, it appears to have a low-passed effect. Why might this be?


title: Conjugate prior for a binomial-like distribution

text: Every week, my $m-1$ friends and I enter a pub quiz which has $n$ points available. In any given week only some of us are there. Record the presence/absence of member $i$ in week $t$ in the matrix $x_{ti}$. I'd like to build a predictive model for the number of points we score in any given week, depending on who's there. A simple model is to assume that each question is equally hard, and that player $i$ has probability $p_i$ of knowing any question. Then the probability that we get a given question correct is $$q({\bf x},{\bf p}) = 1 - \prod_{i=1}^m(1-p_i)^{x_i}$$ Therefore, the probability that we score $k$ points is $$P({\rm Score} = k|{\bf x},{\bf p}) = {n\choose k}q^k(1-q)^{n-k}$$ and the log-likelihood for a given ${\bf p}$ and ${\bf x}$ is $$L = \log {n\choose k} + k\log q + (n-k)\log(1-q)$$ I can write a numerical routine which maximises this, to find the maximum- likelihood estimator for $\bf p$. However, the resulting estimate badly overfits the data (as can be detected by cross validation). A good solution seems to be to introduce a penalty (regularization) term into $L$, which penalizes large or small probabilities. As I understand it, this is equivalent to having a prior on $\bf p$. However, I don't know what the form of this prior should be. Two simple choices for the penalty term are: $$\| {\bf p} - 1/2\|^2$$ and $$-\sum_i\log \left( \frac{p_i}{1-p_i}\right)$$ but these are very ad-hoc. I'd be interested to know what an appropriate conjugate prior is (assuming one even exists). Any hints?


title: R model.matrix and makeContrast. Understanding model and possible contrast

text: I have measurements from 12 mice, grouped in two conditions. I each mouse I have measurements from 4 tissues. The design is not balanced, 5 mice in condition1 and 7 in condition2. After reading the bioconductor edgeR manual I have set up the following model:               design <- model.matrix(~ Condition + Tissue + Condition:Mouse + Condition:Tissue)      I have then manually removed terms with no observations.               design <- design[, -which(colnames(design) %in% c("Condition1:Mouse6", "Condition1:Mouse7"))]      I can fit the model, but when I set up contrasts I get confused. I want to get the main effect of being in condition2, the main effect of each tissue and the interaction effect for each tissue.               colnames(design)      [1] "(Intercept)"  "Condition2"    "Tissue2"                    [4] "Tissue3"  "Tissue4"   "Condition1:Mouse2"                  [7] "Condition2:Mouse2"    "Condition1:Mouse3" "Condition2:Mouse3"               [10] "Condition1:Mouse4"    "Condition2:Mouse4" "Condition1:Mouse5"                  [13] "Condition2:Mouse5"    "Condition2:Mouse6" "Condition2:Mouse7"               [16] "Condition2:Tissue2"   "Condition2:Tissue3"  "Condition2:Tissue4"       To get the main effect of condition I use the contrast               c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      Comparing e.g. Tissue 1 and 2 or 2 and 3 I do like this               c(0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)     c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      But say I want to get the difference between condition 1 and 2 in tissue 1. Given how the model matrix is built tissue 1 is not present. Can I do a comparison like:               c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,-1,-1)      And how about an anova like test for changes across all the tissue:condition interaction terms. The contrast               c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1)      Would compare only tissue2, 3 and 4, and with what? Many questions here. Quick recap. 1\. Is it OK to remove empty terms in an unbalanced design? 2\. How do I find the effect of the level in the tissues category that is missing in the model matrix. 3\. How do I find changes happening across all levels of the interaction terms (related to the missing tissue1 problem). A link to a page with a thorough R model matrix/linear model guide would also be appreciated.


title: Supervised learning algorithm that can be easily retrained with new data

text: I have a web crawler and i want to be able to differentiate a specific class of website (social networks), from others. My problem is that my starting classified data is really small. What I wan't to do is:   * initially train my algorithm on this data;   * receive feedback from user (correct, wrong) on classified websites and use this to increase my data set and retrain my algorithm. Is this feasible, or is my new data biased by having been previously already selected as a social network website by my algorithm? Which supervised learning algorithm would be ideal for being retrained after acquiring new data? Consider that the features of the website I was considering to give the algorithm for training was something like, number of specific html tags (like a [anchor tag], or img) and number of occurrence of some specific words. (If you have suggestions of better features I would be glad to hear also)


title: Anderson Darling exponential distribution

text: I need a goodness of fit test for the exponential distribution. I understand that Kolmogorov-Smirnov is not generally regarded as very powerful and that Anderson-Darling is regarded as superior. However I have two problems.   * None of the literature I can find specifically talks about goodness of fit for exponential distributions. Anderson-Darling is almost always discussed wrt the normal distribution.   * Is there an R (or python) package to do an Anderson-Darling goodness of fit for the exponential distribution? I can find them for other distributions.   * Is there a better test (which exists in R or python)?


title: Transformation for zeros

text: I have been working on seed germination experiments with some dormancy breaking treatments and only for some treatments only I got germination. I need to analyze the parameters like days for initial germination and days for final germination. since I got germination for only few treatments and for others no germination and I need to provide analyzed data regarding initial and final days of germination which is not there for some treatments since no seeds germinated in those treatments. so how to analyze these type of data? Is there any transformation needed?


title: Compare conservation rates of a variable in two-factor experiment

text: This is probably a very simple question for anyone using R, but I am fairly new to it (and stats in general) and I can't figure this simple computation out. I have 3 varieties of plants (A, B, C) tested with 2 different nutrient supplies, and the weights of the plants are measured. The data is balanced. (Here is an example of dataset: http://dl.dropbox.com/u/1755762/example.csv) I can compare the means between genotypes and treatments, with for example `bargraph.CI` from the `sciplot` package, which I like because of the automatic SE bars:               bargraph.CI(Genotype, Weight, Treatment, legend=T, ylab="Dry weight (g)", xlab="Genotype")      which would produce something like that (just to have an idea at what I am looking at): ![enter image description here](http://i.stack.imgur.com/IUb2o.png) I would then complement this with the classes from a Tukey HSD post-hoc test with `HSD.test` from the `agricolae` package (if Anova showed a Genotype×Treatment interaction):               int <- interaction(Genotype, Treatment)     HSD.test(aov(Weight~int), "int", group=T)      But **how can I compare the conservation rates between the two treatments for each genotypes?** By this I mean: compare, for each genotype, the ratio "Treatment 1 weight mean over Treatment 2 weight mean" (this would be just one bar per genotype, showing the weight conservation rate between the two treatments). It would be great to get a graph similar to what `bargraph.CI` produces, but I don't know how to determine if the differences between those weight conservation rates are significant: would standard error bars be relevant here?


title: A problem in probability theory

text: I was trying to solve a problem and I got stuck at the penultimate step (I think). I could show that Var(X) = Var(Y) = Cov(X,Y), where X and Y are random variables with finite means and variances. Based on above statement, can I say that X and Y are the same? or P(X=Y) = 1?


title: Trying to Calculate Basic Lunar Probabilities

text: I'm trying to do what I thought was a basic problem but doesn't seem to be working out properly. I'm looking at disproving a claim that tides cause earthquakes, where the main mechanism that people claim is that it happens during a new or full moon, and especially during a perigee moon (when it's closest to Earth). People like to give ±1 week windows which is basically half a lunar cycle, so it seems basic that the probability an earthquake would happen, given random chance, would be 50% that it would happen within a week of, e.g., a full moon. Similarly, since perigee/apogee happens roughly on the same timescale, it seems as though it would be 50% chance that an earthquake would happen within 1 week of a perigee moon. Put the two together and you have a 25% chance that an earthquake would happen within both a perigee and full moon, correct? (Important note, _updated_ : Actual period between new/full moons averages 29.52 days, actual time between perigee moons averages 27.56 days. **However** , apogee and perigee are NOT normally distributed. Apogee happens a mode of 27.78 days (mean is 27.55±0.27), while perigee is much more asymmmetric, having a mode of 28.4 days at the peak of approximately a Lorentzian. Half the max of the Lorentzian is ±0.08 days. But, the mean is 27.56 with a standard deviation of 1.12; the range is 24.6-28.6 days. I'm thinking this could throw off the modeling?) Assuming that's correct, I seem to be running into issues when trying to figure out the probability that an earthquake would happen by chance within ±X days of BOTH a perigee and full moon. I thought the equation would be simply (2*X/(# days in lunar month, 29.5))*(2*X/(# days between perigee, 27.5)). However, when I do a monte carlo simulation with 500,000 randomly chosen dates within the time period of 1933 to 2012 (just happens to be when I have earthquake data), the fractions do not line up. For example, the simulation shows that 14.4% should be within 5 days of both a full and perigee moon, but my above math says it's only 12.3%. I have checked the results of my simulation against days from perigee, apogee, new moon, and full moon times. As expected from a random distribution with a large $N$, the number of times that the simulated earthquake is a given time period away from maximum perigee/apogee/new/full moon is even. Except for perigee, where I see a fall-off for $>|±12.5|$ days from when it's closest to perigee. I'm thinking this has to do with the non-Gaussian distribution of perigee times? And could that account for the 2% difference at the 5-day example? Is the best way to approach this, because these perigee times are a bit crazy, to simply go with the Monte Carlo results? P.S. This has been updated to better reflect a correction I made in my data. When I initially posted this, I had some incorrect full/new moon dates in my table that were throwing some results off.


title: Is the fixed effect in a fixed effect model a random variable or not?

text: From Wikipedia for fixed effect model: > Consider the linear unobserved effects model for N observations and T time > periods: $$ y_{it} = X_{it}\mathbf{\beta}+\alpha_{i}+u_{it} \quad t=1,..,T, > \quad i=1,...,N $$ where $y_{it}$ is the dependent variable observed for > individual $i$ at time $t$, $X_{it}$ is the time-variant $1\times k$ > regressor matrix, $\alpha_{i}$ is the unobserved time-invariant individual > effect and $u_{it}$ is the error term. Unlike $X_{it}$, $\alpha_{i}$ cannot > be observed by the econometrician. Common examples for time-invariant > effects $\alpha_{i}$ are innate ability for individuals or historical and > institutional factors for countries. > > Unlike the Random effects (RE) model where the unobserved $\alpha_{i}$ is > independent of $x_{it}$ for all $t=1,...,T$, the FE model **allows > $\alpha_{i}$ to be correlated with the regressor matrix $x_{it}$**. Strict > exogeneity, however, is still required. > > Since $\alpha_{i}$ is not observable, it cannot be directly controlled for. > The FE model eliminates $\alpha_{i}$ by demeaning the variables using the > within transformation: $$ > y_{it}-\overline{y_{i}}=\left(X_{it}-\overline{X_{i}}\right) \beta+ \left( > \alpha_{i} - \overline{\alpha_{i}} \right ) + \left( > u_{it}-\overline{u_{i}}\right) = \ddot{y_{it}}=\ddot{X_{it}} > \beta+\ddot{u_{it}} $$ where > $\overline{X_{i}}=\frac{1}{T}\sum\limits_{t=1}^{T}X_{it}$ and > $\overline{u_{i}}=\frac{1}{T}\sum\limits_{t=1}^{T}u_{it}$. Since > **$\alpha_{i}$ is constant** , $\overline{\alpha_{i}}=\alpha_{i} $ and hence > the effect is eliminated. The FE estimator $\hat{\beta}_{FE}$ is then > obtained by an OLS regression of $\ddot{y}$ on $\ddot{X}$. On one hand, it says "the FE model allows $\alpha_{i}$ to be correlated with the regressor matrix $x_{it}$", which I understand to be that both $\alpha_{i}$ and $x_{it}$ are random variables. On the other hand, it says "$\alpha_{i}$ is constant", and my previous impression is that $\alpha_{i}$ is nonrandom in fixed effect model, and random in random effect model. So I was wondering how to correctly understand "fixed effect"? Thanks and regards!


title: What test should be used for detecting team imbalances in a game?

text: I'm playing the iPhone/Steam game Hero Academy. It's a bit like Chess, except there's some randomness, and at the beginning of a game players choose what team they'll use, each with different strengths and weaknesses. Both players can pick the same team. About 400 of us play in an unofficial league, where we record our games online for Elo rankings. I'm not sure how many are active players, though. I have the win/loss record for every matchup between the four teams. Draws are very rare. The rarest matchup has 143 games; the most-common matchup has 260. There are 3,234 games in the database. **What is the appropriate test for determining whether the teams are of equal strength?** Let's assume that rock-paper-scissors imbalances are undesirable. That is, if Team A beats Team B which beats Team C which beats Team A, we have an imbalanced game. Every matchup should be fair. A complication is that many players have preferred teams, so that if Team A is very popular among top players, it'll appear stronger. I suspect we'll have to ignore this effect for now. Another complication is that there's a slight advantage to going first, despite our efforts to ameliorate this with a handicap. Let's ignore this effect as well.


title: SVM, Overfitting, curse of dimensionality

text: My dataset is small (120 samples), however the number of features are large varies from (1000-200,000). Although I'm doing feature selection to pick a subset of features, it might still overfit. My first question is, how does SVM handle overfitting, if at all. Secondly, as I study more about overfitting in case of classification, I came to the conclusion that even datasets with small number of features can overfit. If we do not have features correlated to the class label, overfitting takes place anyways. So I'm now wondering what's the point of automatic classification if we cannot find the right features for a class label. In case of document classification, this would mean manually crafting a thesaurus of words that relate to the labels, which is very time consuming. I guess what I'm trying to say is, without hand-picking the right features it is very difficult to build a generalized model ? Also, if the experimental results don't show that the results have low/no overfitting it becomes meaningless. Is there a way to measure it ?


title: How to use test set data if model has been built using a training set transformed with PCA?

text: I didn't know this forum existed so I asked this on Stack Overflow but maybe it's more suitable here... I have a question regarding dimensionality reduction using PCA. If I have a train set (train.arff, 10 attributes) I perform a PCA and I save my data with respect to the new transformed variables (say I choose the two first attributes, combination of the original ones, that collect most of the variance), and call this transformed trainset "trainset-afterPCA.arff". Now I train a model using this file (which only has 2 attributes), and save it. If now I have a new dataset, constructed with the original 10 attributes, and I want to use the model I built before to classify this new data, how do I have to proceed? If I just try to test on this new dataset, train and test aren't compatible, right? If I ran PCA on the test set, the resulting new attributes won't be the same as the ones obtained in the training set. What should I do? * * * Update: Just to show (part) of the weka output:               eigenvalue  proportion  cumulative     2.31715   0.28964     0.28964   0.512ent-0.472Threshold+0.422impRes-0.335pssm-mut+0.28 pssm-wt...     1.72263   0.21533     0.50497   0.593pssm-mut+0.501pssm-wt+0.41 Threshold+0.403hyd+0.161sub...     1.31987   0.16498     0.66996   0.698vdw+0.628sub+0.219Threshold-0.168hyd+0.154ent...     0.88362   0.11045     0.78041   0.53impRes-0.51pssm-wt-0.478ent+0.346hyd+0.33 subs-score...     0.8404    0.10505     0.88546   0.605hyd-0.552impRes-0.319pssm-wt-0.26pssm-mut+0.235ent...     0.56935   0.07117     0.95663   -0.656vdw+0.531sub-0.449hy+0.207Threshold-0.145pssm-mut...          Eigenvectors     V1       V2       V3     V4      V5      V6      -0.4716  0.4104  0.219  -0.0231  0.215   0.2071 Threshold     -0.153  -0.1263  0.6977  0.0049 -0.1865 -0.6556 vdw      0.2465  0.4028 -0.1679  0.346   0.6055 -0.4486 hyd      0.2511  0.1609  0.6277  0.3299  0.1513  0.5306 sub      0.2799  0.5007  0.0529 -0.5097 -0.3189  0.061  pssm-wt     -0.335   0.593  -0.1004  0.0423 -0.2602 -0.145  pssm-mut      0.5118  0.0616  0.1544 -0.4783  0.2354 -0.1246 ent      0.4217  0.1459 -0.0799  0.5297 -0.5521 -0.0645 impRes      So, based on your answer, if now I choose V1 and V2 to represent my data, I have to use V1 and V2 and calculate the new attributes for the test set, previous to upload the test set to the model..                V1 --> new att 1 = 0.512ent-0.472Threshold+0.422impRes-0.335pssm-mut+0.28pssm-wt...      V2 --> new att 2 = 0.593pssm-mut+0.501pssm-wt+0.41Threshold+0.403hyd+0.161sub...


title: Inference from conditional observations

text: Let $(x_1, \ldots, x_n)$ be an i.i.d. random sampling from a conditional normal distribution ${\cal N}(\mu,\sigma^2)$ distribution given some event $A$ possibly parameter-dependent: for instance when we sample ${\cal N}(\mu,\sigma^2)$ but we only retain $x_i$ when $x_i \in [\mu-3\sigma, \mu+3\sigma]$. Are there some known methods to perform statistcal inference on the parameters?


title: Out of sample predictions for logistic regression models in R

text: I have a problem of the form "what is the probability that a user will 'like' a certain movie?" For a bunch of users, I know the movies each has watched historically, and the movies each has liked. Additionally, for each movie I know the name of the director. I calibrated a logistic regression for each user of the form: `glm(liked_by_user_1 ~ liked_by_user_2 + ... + liked_by_user_k + factor(director), family=binomial, data = subset(MovieWatchings, user_id == 1))` But my problem is: say that in the past, user 1 has watched movies from directors `D1` through `DM`, but next month `U1` watches a movie directed by `DN`? In that case the R `predict()` function will give an error, because the glm model for user 1 doesn't have an estimated parameter for the case of `director = DN`. But I must know something about `U1's` probability of liking the new movie, because I still know which other users have seen and liked this movie, and that has some predictive power. How can I set up my model so that I can take into account other users' liking behavior, AND user 1's director preferences, but still have sensible predictions when user 1 sees his first movie from a new director? Is logistic regression even the right type of model for this case?


title: How to find best split in decision trees using label vectors?

text: In single label learning, we have a single label that we need to predict. Most training tasks are single label, like the most common example where we need to predict if a customer is going to be good or bad. In multi label learning there are more than one labels. For example if we want to tag a photo, we put more than one tags. The labels now are described by a vector and not by single values like in single label learning. I am trying to build a decision tree that finds best splits based on variance. Me decision tree tries to maximize the following formula:                  Var(D)*|D| - Sum(Var(Di)*|Di|)      D is the original node and Di are the splits produced by choosing an attribute (by Di, i mean the node that is produced by choosing an attribute and its i-th value). My problem is this. Imagine that i have a matrix for storing each node's label vectors, where in rows are the examples and in columns the values for every label. Lets say that this matrix is mXn. After finding the variance for this matrix, i have an 1Xn vector, where in each column, i have the variance for every label. Now i need to find the the difference between 2 vectors, where each vector is multiplied by a number. My problem is how to decide which is the best split, when the outcome of my formula is a vector? I mean that if i had single labels, i would have a single number as an outcome and i would simply choose the biggest one. Now that i have vectors, how can i decide which is the best split? For example lets say that the outcome from my formula for 2 different attributes, are these 2 vectors:                1,2,3,4,5    (attribute i)      2,3,4,5,6    (attribute j)      How do i choose between attribute i and j?


title: Seasonally adjusted month-to-month growth with underlying weekly seasonality

text: As a side hobby, I have been exploring forecasting time series (in particular, using R). For my data, I have the number of visits per day, for every day going back almost 4 years. In this data there are some distinct patterns:   1. Monday-Fri has a lot of visits (highest on Mon/Tue), but drastically less on Sat-Sun.   2. Certain times of the year drop (i.e. many less visits around U.S. Holidays, summers show less growth)   3. Significant growth year-to-year It would be nice to be able to forecast an upcoming year with this data, and also use it to have seasonally adjusted month-to-month growth. The main thing that throws me off with a monthly view is:   * Certain months will have more Mon/Tue than other months (and that isn't consistent over years either). Therefore a month that happens to more weekdays needs to be adjusted accordingly. Exploring weeks also seems difficult since the week numbering systems change from 52-53 depending on the year, and it seems `ts` doesn't handle that. I'm pondering taking an average for the weekdays of the month, but the resulting unit is a bit strange (Growth in Avg. Weekday Visits) and that would be dropping data which is valid. I feel this sort of data would be common in time series, (say for example electricity usage in office building might be something like this), anyone have any advice on how to model it, in particular, in R? The data I am working with is pretty straight forward, it starts like:                           [,1]     2008-10-05 17607     2008-10-06 36368     2008-10-07 40250     2008-10-08 39631     2008-10-09 40870     2008-10-10 35706     2008-10-11 18245     2008-10-12 23528     2008-10-13 48077     2008-10-14 48500     2008-10-15 49017     2008-10-16 50733     2008-10-17 46909     2008-10-18 22467      and continues like this up to the present, with an overall trend of growth, some dips around US holiday weeks, and growth generally slowing during the summer.


title: How to get started with rating and ranking based on pairwise competition data?

text: I'm interesting in learning about how to rate and rank individuals in a group who only interact/compete in a pairwise fashion (i.e., systems like the ELO rating system for chess).   * Are there are any go-to methods or more accurate and advanced methods out there?    * Are there any R packages that make implementation easy?    * Are there methods that can use auxiliary information as well as the outcome of a match/game?    * Are there methods that can better use the information of winning margin as opposed to dichotomous win/lose?   * What should I be looking for in the literature?


title: Random number-Set.seed(N) in R

text: I realize that one uses set.seed() in R for pseudo-random number generation. I also realize that using the same number, like set.seed(123) insures you can reproduce results. But what I don't get is what do the values themselves mean. I am playing with several functions, and some use set.seed(1) or set.seed(300) or set.seed(12345). What does that number mean (if anything)- and when should I use a different one. Example, in a book I am working through- they use set.seed(12345) when creating a training set for decision trees. Then in another chapter, they are using set.seed(300) for creating a Random Forest. Just don't get the number. Thanks


title: How to add non-linear trend line to a scatter plot in R?

text: I have a scatter plot. How can I add non-linear trend line?


title: Understanding the linear mixed effects model equation and fitting a random effects model with weights in R

text: Could anybody please shed some light on me? I have $N$ individuals, divided into $K$ groups with $N>K$. Some groups have 1 individuals only, some other groups have more individuals. Each individual has 4 features (4 variables): $F_0, F_1, F_2,$ and $F_3$, where $F_0$ is simply the group indicator. Therefore the data matrix is of size $N \times 4$. There is also a weight vector $w$ of length $N$, which gives the weights for each of the $N$ individuals in the regression. May I ask if the following model is a random-intercept model?   1. There is a common beta for all N individuals.   2. Each group has a different within group regression line (same slope but different intercepts).   3. The regression line within each group crosses the "cloud" consisting of the group members. And the individual residuals scatter around the regression line, within each group. This sounds like a "random-intercept" model to me. However, how do I explicitly write out the equation? $$ {\boldsymbol y} = {\bf X} \beta + {\bf Z} b + \varepsilon $$ More specifically, with the three feature variables $F_1, F_2$ and $F_3$ and the group indicator variable $F_0$. I am having difficulty writing out $X$ and $Z$ explicitly. Moreover, my $F_2$ is a factor variable. Could anybody please show me how the $X$ and $Z$ matrices look like explicitly? And what do the "$b$"'s represent here? And how do I set up the weights in LME in `R`? If I would like to have "group"-weights and "individual"-weights, how shall I do it? Thank you! [ps1 @ Macro] Thanks a lot Macro for your very comprehensive answer. I am continuing digesting your writings and I will consult with you and your expertise with more questions. My first question is: you think the condition 2 and condition 3 conflict each other. I don't understand why. Here is my thinking: even though for a certain group, the overall beta is a very crappy one for that group. By adjusting the intercept for that group, I at least still can get the regression line passing the cloud, am I right? More specifically, in a random intercept model (all group sharing the same beta), will the following situation arise? All data points in a group reside on one side of the regression line of that group? * * * In a random intercept model, for each group, which point does the regression line pass thru? In the general OLS, we knew that it's the point (xbar, ybar) that is passed thru by the regression line. But for each group in a random intercept model, which point is that "central" point? * * * Moreover, in R, what's a convenient way to visualize what's happening within a certain group in a mixed model? Thanks a lot for your help!


title: Dropped cases from matched studies

text: We have cohort data and a rare exposure which we are matching to controls in a large epidemiologic dataset. The matching variable is a deidentified neighborhood indicator (cluster) which guarantees that individuals within clusters are highly correlated. This is a social epidemiology type of analyses looking at health related outcomes and socio-behavioral "determinants" as possible risk factors with cross-sectional data (our analytic variables of interest were only collected at baseline). After controlling for individual level confounders, we were left with the issue of neighborhood level confounders which are vast and of extreme importance. Neighborhood level information was not made available to us, so matching based on neighborhood seems like an adequate way to control for its confounding effects. When the matching was performed, we lost 10% of cases due to unmatched pairing (small numbers or lack of controls). My question is:   1. does exclusion of unmatched cases indicate a possibility for residual confounding due to unmeasured neighborhood level confounders (despite all individuals being adequately matched in the final analysis)   2. is looking at a fixed effects model a feasible sensitivity analysis for assessing whether matching was reasonable and correct?


title: pvals.fnc does not work with random intercept and slope. How to get p-values?

text: I am doing LMM in lme4, and have installed languageR to run pvas.fnc(), but it does not work with random intercept and slope. This is the message:               Error in pvals.fnc(lmer1) :        MCMC sampling is not yet implemented in lme4_0.999375       for models with random correlation parameters      For model with only random intercept, this function works fine. Then what should I do to get the p-values?


title: What are my chances of winning?

text: Let's assume a hypothetical contest. You can join the contest by submitting a ticket, but multiple tickets can also be submitted by the same person. (For the purpose of this question, let's say there are 5000 tickets entered and 300 applicants). Twenty applicants can win, no person can be selected twice. What are the chances of one person winning with n tickets submitted?


title: Compare LMM GLMM (generalised linear mixed model, negative binomial) by numerical measure (AIC BIC, cross validation, R² squared) for model validation

text: How to compare results of generalized linear mixed model (GLMM, negative binomial) with a log transformed linear mixed model (multilevel, hierarchical) . I have a data set (counts), which is nested. I fitted 2 mixed models: (1) linear mixed model with a log transformed response variable and (2) negative binomial mixed model. My question is: How to compare these two models? Which one is the "better" model? How much "better" is the "better" model? I would need a comparable measure. As I have no zeros and a large mean, log transformation should approximate the normal distribution quite well. I know, there is some rejection of transforming count data in literature (see eg. Hara & Kotze. 2010. Do not log- transform count data). However, results from these simulated data may not be transferable to real world applications (as stated by the authors them self). Still, I wonder which one is the "better" model and to which extend it is better. My validation plots do not allow to identify a preferable model. I know, that this is a controversial issue and I'm not expecting a simple black or white answer. However, are there some pragmatic approaches? The data set:   * about 220 observations (counts),    * 12 predictors in the full model,   * mean is relatively large and   * large variance (Poisson does not work due to high overdispersion)   * no zero observations   * 2 random slopes The fitted models (after stepwise variable selection): (1) linear mixed model with log transformed response                 LMM <- lme (log(Y) ~ a + b  + e + f, random = ~ - 1 +  e + f | location        , method = "REML", data=my_data)                                          (2) negative binomial mixed model                 GLMM  <- glmmPQL(Y ~ a + c + d + e + f, random = ~ - 1 +  e + f | location       ,data=my_data)      (A) Some sort of R² (R squared) could be calculated LMM and GLMM (marginal & conditional). See for eg. http://glmm.wikidot.com/faq and here (Nakagawa & Schielzeth (2013) A general and simple method for obtaining R 2 from generalized linear mixed-effects models). However, I could not apply it on a negative binomial mixed model. (B)Information criteria (AIC, BIC etc.) can be calculated. They are also debated for mixed models (http://glmm.wikidot.com/faq ), in particular how degrees of freedom are to be calculated. I ran into the following difficulties:   * using glmmPQL from the MASS package for fitting the GLMM does not produce AIC scores (glmmPQL uses quasi-likelihood and no log-likelihood (hence no AIC nor BIC, which are based on the log-likelihood)   * other packages for fitting GLMM with negative binomial distribution such as glmmADMB or glmer.nb function from the lme4 package produces AIC scores, but on a very different scale as compared to the lme function from the nlme package for fitting the LMM. Furthermore, I had difficulties fitting my GLMM using glmmADMB or glmer.nb (convergence and some strange error messages)  (C) Cross validation:   * I did not manage to run cross validation for my GLMM (negative biniomial) using glmmPQL Does anyone know how to estimate those measures to evaluate the 2 models? Or are there any other comparable measures? Any hints are very welcome! R-code is appreciated as well. Thanks Best Philipp


title: When training a RBM how should the data be normalised?

text: A feed forward neural network trains best if the data is normalised so that each input has a mean of 0 and a standard deviation of 1. Is this true for a Restricted Boltzmann machine as well? (My experiments are saying, no it isn't.)


title: How do I test if two (non-normal) distributions differ?

text: I have read about Student's t-test but it appears to work when we can assume that the original distributions are normally distributed. In my case, they are definitely not. Also, if I have 13 distributions, do i need to do `13^2` tests? ![Here is a sample of two distrbutions.  There are 13 distributions.](http://i.stack.imgur.com/k6Smy.png)


title: How to get p-value for fixed-effects by model comparsion?

text: I am doing linear mixed model using lme4. According to Winter (2013, http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf), as the new version of R does not give p-values due to inconclusion of degree of freedom, p-values of mixed models can be derived from model comparison. From the examples on page 12, he suggested to construct the null model:               politeness.null=lmer (frequency ~ gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)      Then add the fixed effect that we are interested in:               politeness.null=lmer (frequency ~ attitude + gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)      And then the p-values of attitude can be given from               anova(politeness.null,politeness.model)      However, in my case, I have 3-way interaction: color*sex*food, and when I run the model I have 17 layers of fixed effects, such as white (compare to red), blue (compared to red), white:male(compared to female), etc. Then my question is how I can get p-values for all these fixed effects? I am not sure if I should have some fixed effects that I am not interested in first:               lmer1=lmer (duration ~ action + (1|subject) + (1|repetition), data=data.frame, REML=FALSE)       Then add:               lmer1=lmer (duration ~ action + color + (1|subject) + (1|repetition), data=data.frame, REML=FALSE)        Then add:               lmer1=lmer (duration ~ action + color + sex + (1|subject) + (1|repetition), data=data.frame, REML=FALSE)       or               lmer1=lmer (duration ~ action + color*sex + (1|subject) + (1|repetition), data=data.frame, REML=FALSE)      or               lmer1=lmer (duration ~ action + sex + (1|subject) + (1|repetition), data=data.frame, REML=FALSE)      Anyone can help me about this? or is there any other ways to get p-values easier than model comparison?


title: How is "Orthogonal distance" computed?

text: I was reading the vignette of the R package `chemometrics` (link). In the second paragraph (right below the first equation) of Page 12, the author writes: > the OD (Orthogonal Distance) is calculated in the original space as the > orthogonal distance of an object to the PCA subspace or, in other words, the > distance between the object and its orthogonal projection on the PCA > subspace. Can someone elaborate upon how exactly the orthogonal distance is calculated (some R code as an illustrative example would be greatly appreciated, if possible!)


title: Interpreting random slopes equal to 0 using lmer in R

text: I have just started working with mixed models and the lme4 package and was after some advice interpreting some results. I have a data set looking at the change in nest height of birds (NAP) over time (Year). I have two random effects that are relevant to my data; the area where the birds nest (Area) and the identity of the birds (MaleID). I started by fitting a simple mixed model with random slopes for MaleID, as this is most relevant to my question.               ELEV.LME<-lmer(NAP~Year+(0+Year|MaleID),data=ELEV)      This model works fine and shows some variance in the slopes of MaleID. However, I also know that Area and possibly random intercepts of MaleID are also important. I fitted a new model including area;               ELEV.LME2<-lmer(NAP~Year+(0+Year|MaleID)+(1|Area),data=ELEV)      but now the variance for random slopes in MaleID is 0. The same issue occurs if I include random intercepts for MaleID;               ELEV.LME4<-lmer(NAP~Year+(1|MaleID)+(0+Year|MaleID),data=ELEV)      I originally thought this may be caused by having cases of MaleID that only occur once in the data set, but I received the same 0 value when these were removed. Is this 0 value caused by my coding or a structural issue with my data? Or does this result tell me that variation in the slopes of MaleID explains none of the variance in NAP over Years? I wasn't expecting a large variance in random slopes, but 0 seems quiet strange, especially as I received the same result after I change around the dataset. Any help or advice would be greatly appreciated. Thank you! I have included a dput() output from a subset of my data below. The real data set is thousands of records, and has many more levels of Area, but I encounter the same result with this reduced data set.               structure(list(MaleID = c(1003L, 1010L, 1113L, 1113L, 1113L,      1113L, 1113L, 1113L, 1113L, 1113L, 1113L, 1113L, 1113L, 1113L,      1167L, 1167L, 1167L, 1168L, 1168L, 1168L, 1200L, 1207L, 1207L,      1216L, 1216L, 1216L, 1216L, 1216L, 1216L, 1216L, 1223L, 1223L,      1240L, 1240L, 1240L, 1240L, 1240L, 1240L, 1240L, 1257L, 1257L,      1299L, 1300L, 1300L, 1301L, 1302L, 1302L, 1310L, 1312L, 1327L,      1327L, 1327L, 1330L, 1347L, 1347L, 1347L, 1347L, 1347L, 1347L,      1347L, 1347L, 1347L, 1379L, 1390L, 1393L, 1393L, 1393L, 1393L,      1393L, 1393L, 1393L, 1393L, 1404L, 1404L, 1420L, 1420L, 1420L,      1420L, 1420L, 1420L, 1425L, 1440L, 1440L, 1449L, 1449L, 1449L,      1449L, 1449L, 1449L, 1449L, 1449L, 1449L, 1452L, 1452L, 1454L,      1454L, 1465L, 1465L, 1465L, 1465L, 1466L, 1488L), Year = c(1995L,      1995L, 1995L, 1996L, 2008L, 2010L, 2010L, 2011L, 2011L, 2012L,      2012L, 2012L, 2013L, 2013L, 1995L, 1996L, 1996L, 2008L, 2010L,      2011L, 1995L, 1995L, 1996L, 2008L, 2010L, 2011L, 2011L, 2012L,      2012L, 2013L, 1995L, 1995L, 2008L, 2008L, 2010L, 2011L, 2011L,      2012L, 2013L, 1996L, 2008L, 1995L, 1995L, 1996L, 1996L, 1996L,      1996L, 1996L, 1996L, 1996L, 2008L, 2010L, 1996L, 2008L, 2008L,      2010L, 2011L, 2011L, 2012L, 2012L, 2012L, 2013L, 1996L, 1995L,      1995L, 1995L, 1996L, 2008L, 2008L, 2010L, 2011L, 2011L, 1995L,      1995L, 1996L, 1996L, 2008L, 2010L, 2010L, 2010L, 1996L, 1995L,      1995L, 1995L, 1996L, 2008L, 2010L, 2010L, 2011L, 2011L, 2011L,      2012L, 1995L, 1996L, 1996L, 1996L, 1995L, 1995L, 1995L, 1995L,      1996L, 1995L), Area = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L), .Label = c("A",      "B"), class = "factor"), NAP = c(158L, 144L, 144L, 150L, 149L,      155L, 143L, 140L, 147L, 140L, 132L, 143L, 137L, 141L, 147L, 155L,      138L, 140L, 190L, 163L, 144L, 139L, 138L, 144L, 147L, 142L, 135L,      136L, 141L, 144L, 153L, 149L, 135L, 143L, 209L, 144L, 220L, 138L,      185L, 135L, 136L, 189L, 182L, 145L, 189L, 138L, 184L, 149L, 150L,      157L, 178L, 187L, 190L, 147L, 145L, 149L, 145L, 166L, 140L, 138L,      144L, 143L, 158L, 138L, 140L, 135L, 142L, 134L, 152L, 153L, 142L,      198L, 169L, 142L, 151L, 160L, 160L, 163L, 153L, 153L, 141L, 136L,      148L, 142L, 175L, 153L, 144L, 162L, 168L, 169L, 142L, 154L, 139L,      129L, 161L, 166L, 162L, 134L, 159L, 152L, 138L, 152L)), .Names = c("MaleID",      "Year", "Area", "NAP"), class = "data.frame", row.names = c(NA,      -102L))


title: Comparing two Naives Bayes Classifiers (using different features of the same data) and using incremental k-fold cross-validation

text: > **Possible Duplicate:**   >  Comparing two classifier accuracy results for statistical significance with > t-test I coded two Naives Bayes Classifiers (using different features of the same data) and used incremental k-fold cross-validation. As output I have computed (for each of the two NBC): For each k training set sizes: Average accuracy / Standard Error Upon observing this data I formulate the hypothesis that one of the two has better performance (in terms of accuracy). How can I assess the significance of my results in terms of whether I should reject the Null hypothesis in favor of my hypothesis or not? Thanks for any help :)


title: Path analysis, sample sizes, and alternative analysis

text: I am examining how English ivy affects the occurrence of a salamander species under cover objects (e.g., logs). Soil moisture is assumed to be the major factor that affect their occurrence. My hypothesized pathway: The presence/absence of salamanders under cover objects is either a direct consequence of changes in ivy-induced abioitc environment (i.e., drier soil) or an indirect result of changes in prey community that resulted from altered abiotic factors. But, there are multiple factors, other than English ivy, that affect soil moisture. ![enter image description here](http://i.stack.imgur.com/k65Ag.jpg) My questions are:   1. I think that a path analysis is most suitable for testing my causal mechanisms. But, given a small sample size (n = 71), is a path analysis appropriate?   2. Another potential problem for a path analysis is that the effects of English ivy on soil moisture seem to depend on the other factors (e.g., the number of overstory trees), as shown below. Are there any way to account for such patterns in a path analysis? ![The relationship between soil moisture and English ivy cover on cover objects \("the number of overstory trees" for the left graph\) for different levels of the surrounding overstory trees \("English ivy cover on cover objects" for the left graph](http://i.stack.imgur.com/ArgZm.jpg)   3. Are there any other analyses suitable for testing my hypothesized relationships? I am considering multiple (linear and logistic) regressions, but again my sample size is small **AND** regressions do not reflect my hypothesized causal relationships accurately. I am using R, so any recommended code would be greatly helpful (I am a relatively new R user, though).


title: STATISTICA 12 for SEM data analysis

text: I'm a graduate student pursuing Ph.D. in Information Systems. My dissertation research involves using structural equation modeling (SEM-PLS) as a main data analysis method. After comparing various software packages, I came to conclusion that the optimal solution for me would be to rent STATISTICA 12 (Ultimate Academic Bundle) in order to perform data analysis (basic statistical analysis and SEM) after data collection. I have the following questions:   1. What are advantages and disadvantages of using STATISTICA's 64 bit version versus the 32 bit version?   2. Does STATISTICA 12 have a facility to perform and document pilot and main SEM studies? Thank you in advance!


title: Training one class SVM using LibSVM

text: I hope to use one-class SVM of LIBSVM to train a training samples so as to get a model. Then, I will use the model to predict whether the new test data and the training data is same type or not. In the training process, I have some questions as follows:   * Should the training samples all be positive examples or not?   * Which kernel function can get better result, **linear** kernel or **RBF** kernel?   * What is the effect of nu's values to the model?


title: Method to determine distribution of data

text: Is there a book or paper which mentions the ways how to find the distribution of observed signals. After the mean and variance is found out, how do I determine/find the distribution to which it belongs? There are software packages to do that for you but is there any theoretical framework? I have googled and found tests like chi-square and all but there is no analytical form of calculating mean and other parameters and where it mentions how to find the distribution of data


title: What will be the correct answer, if we modify the "Best statistics question ever"?

text: There is a popular question, called "Best statistics question ever". If you choose an answer to this question at random, what is the chance you will be correct? A) 25% B) 50% C) 60% D) 25% This task is not very difficult, the correct answer is 0%. But if we modify it like this: If you choose an answer to this question at random, what is the chance you will be correct? A) 50% B) 25% C) 60% D) 50% What will be the correct answer? Do we have two correct answers: 25% and 50%, or there is no correct answer, as with this two correct answers the chance of choosing the correct answer is in fact 75% (but we do not have 75% written on the desk)? By the way. Does the answer 0% remains correct answer, the third correct answer in this case?


title: Difficulty with MCMC implementation

text: I could really use some guided help! I'm having difficulty understanding an MCMC implementation in terms of modeling a data set. I'm working on generating parameters from stellar light curves, and was asked to look into implementing an MCMC algorithm. A large chuck on the code is written in Python, so I've been trying to use emcee hammer to generate parameter fits. But going through the code, it's just not "clicking" how the method works. I have a set of data (time vs flux) of two stars orbiting each other such that from our point of view, they eclipse. There are dips in the light curve to signify this. All I'm attempting to do is get the parameters of the system dependent on the characteristics of these dips. In the emcee implementation, there are a few functions that I understand: the posterior function which, I believe, simply generates a data set given the set of parameters. Then there's a prior function which, I assume, is the function given a previous set of parameters. Somehow the algorithm chooses whether or not the jump to the posterior parameter set is to be done? I'm guessing that's what the use of the likelihood function is? To describe whether or not to take the jump? I apologize, I'm quite confused on how this is to be implemented in terms a defined set of data.


title: Exponential family parameter estimation and fitting, references

text: First of all, I want to express my apologies if the question is too broad or wrong, but I am in need of references and I have no idea whom I can ask. If you are interested, the question comes from a model I built, you can see some details here and here. In this model I have: $$f(\mathbb{x}|T,\mu)=\frac{h(\mathbb{x})e^{-\frac{E(\mathbb{x})}{kT}+\mu N(x)}}{\mathcal{Z}(T,\mu)}$$ There, my parameters are $\mu$ and $T$, and $\mathbb{x}=(x_1,\dots,x_M)$ where $x_i\in\\{0,1\\}$ and I have the restriction $\forall i\in\\{1,\dots,M-D+1\\}$ $$\sum_{j=0}^{D-1} x_{i+j} \leq 1$$ This is, $h(\mathbb{x})=0$ if that condition is not held. I have the "small" inconvenience of not knowing $\mathcal{Z}(T,\mu)$, so I used a MCMC (Metropolis-Hastings) method to approximate this function. However I face two problems.   * The first of them regards the simulation and the model and I am on solving it (it depends too much on the initial condition).    * The second is that these parameters are not fully known and I have no idea how can I estimate them. I have been reading about Bayesian inference and I know a bit of estimation theory but I am no expert (furthermore I don't know if not knowing the partition function can affect the result). If any of you were able to give me some clue in the form of a book that I can read, I would be eternally grateful. Thank you very much for your help. Thanks to cardinal's comment, I have realized that I didn't explain one thing. It probably makes all more complex but there it goes: The idea is that $E$ is known in each experiment, actually $E(\mathbf{x}) = \mathbf{E}\cdot\mathbf{x}$. However, $\mathbf{E}$ is not always the same, it represents an external potential for some particles. The "good" thing is that $T$, which accounts for the temperature, never changes whatever $\mathbf{E}$ is, so I thought that I could find a way of estimating it, given the fact that I have an empirical distribution of $x_{i}$ (so, a probability that a particle is in the position $i$) given a certain $\mathbf{E}$. So, in a way, what I have is $$f(\mathbf{x}|T,\mu , \mathbf{E})$$, but I always know $\mathbf{E}$ and I know (can I say this?) that $T,\mu$ are independent of $\mathbf{E}$. I am sorry for not being clear enough before. I am starting to think that nothing of this makes sense...


title: Forced-choice data possible with confirmatory factor analysis?

text: I have a task where I ask participants to view a series of images. After viewing each image they are given several options with which they must a forced-choice as to which category the image falls within. The categories are often quite close to each other conceptually or visually e.g. different ethnicities from a face (e.g. white, asian, black), different facial expressions (e.g. happy, angry, sad), different periods for a piece of architecture (e.g. neoclassical, georgian, gothic. My question is this: does such data lend itself to analysis by exploratory and confirmatory factor analysis (using the sum of "correct" responses for each category), or is there a problem inherent in the somewhat non-independent nature of these measures that makes this analytic approach problematic? Grateful for any help or pointers.


title: 1-step-ahead predictions with dynlm R package

text: I've fit a model with several independent variables, one of which is the lag of the dependent variable, using the dynlm package. Assuming I have 1-step-ahead forecasts for my independent variables, how do I get 1-step-ahead forecasts for my dependent variables? Here is an example:               library(dynlm)          y<-arima.sim(model=list(ar=c(.9)),n=10) #Create AR(1) dependant variable     A<-rnorm(10) #Create independant variables     B<-rnorm(10)     C<-rnorm(10)     y<-y+.5*A+.2*B-.3*C #Add relationship to independant variables      data=cbind(y,A,B,C)          #Fit linear model     model<-dynlm(y~A+B+C+L(y,1),data=data)          #Forecast     A<-c(A,rnorm(1)) #Assume we already have 1-step forecasts for A,B,C     B<-c(B,rnorm(1))     C<-c(C,rnorm(1))     y=window(y,end=end(y)+c(1,0),extend=TRUE)     newdata<-cbind(y,A,B,C)     predict(model,newdata)      And here is an example using the dyn package, which works.               library(dyn)          #Fit linear model     model<-dyn$lm(y~A+B+C+lag(y,-1),data=data)          #Forecast     predict(model,newdata)the dyn packages, which works:


title: Removing duplicated rows from R data frame

text: Coulmn1 Column 2          A    1          A   1          A   2          B   4          B   1          B   1          C   2          C   2      I want to remove the duplicates based on both the columns so that my final output will be               A   1          A   2          B   4          B   1          C   2      order is really not a matter for me Thank in advance


title: When to use Equal-Frequency-Histograms

text: ... instead of e.g. the popular Equal-Width-Histograms. Additional question: What is a good/robust rule of the thumb to calculate the number of bins for equal frequency histograms (like the Freedmann-Diaconis- Rule for equal-width).


title: What should I do to compare different sets of data?

text: I am a beginner in statistics, and I want to learn machine learning :). Therefore, I have gathered some sample data to practice. But, the problem is I want to create a feature (or attribute), which is common for every entity. But, this feature is observed and measured for entity X in Y meters and Z seconds, and for entity W in M meters and L seconds, and there are many observations like described before. Which algorithm should I apply to create a common feature for all entities? BTW, I am using WEKA.


title: Reading in SVM files in R (libsvm)

text: The data files from http://www.csie.ntu.edu.tw/~cjlin/libsvm/ are in 'svm' format. I am trying to read this in to sparse matrix representation in R. Is there an easy/efficient way to do this? Here is what I am doing now: read in file line by line (800,000 lines), for each line separate classes, values, and cols. Store the classes as a list and the features as a .csr sparse matrix (1 row), then rbind the feature row with all previous rows. This is terribly inefficient and basically won't finish (12 minutes for 1000 lines). I think it comes from rbinding the sparse matrices once the number of rows starts to get large. Note: the matrix (800000*48000) is too big to build and then convert to sparse format. Thanks!


title: Comparing groups using SPSS

text: I have a dataset with 2 measures that were tested over 30 time points. I have a two groups of people who have both measures at all 30 time points. How do I compare the groups using SPSS?


title: Monte Carlo integration

text: I am calculating a simple integral $\int^1_{-2} \exp{x^2}(x+1)dx $ with Monte Carlo method using a linear density function $p_\xi (x) = \frac{4}{9} + \frac{2}{9} x $. Let say I have a a sample which is distributed on the interval [0,1]. In order it to be distributed on interval [-2,1] I integrated it, following the algorithm of calculation, $\int\limits_0^{\xi} \frac{4}{9} + \frac{2}{9} x dx= U $ After solving the equation I get $\xi = \sqrt{9 \cdot U + 4} – 2$ . As ξ should be from [-2,1], I transformed a little bit the equation, so my final result was $\xi = \frac{3(\sqrt{9 \cdot U + 4} – 2 )}{\sqrt{13} – 2} – 2 $, where U is random variable from my sample. The rule says that the approximation of integral should be $I \approx \frac{1}{n} \sum \limits^n_{i=1} \exp({\xi_i^2})(\xi_i +1 ) \cdot \frac{1}{p_{\xi} (\xi_i)} $ Where n is a number of random variables in U. I have tried to calculate it with R, but I always get a wrong answer. The true result is -8.02465, but I get -1163.554. I don‘t know where exactly my mistake is, whether it is in my R code or the algorithm. R code:               x1=c()     x1[1]=1     a=1141     c=247     b=a-1     m=1216     for (i in 1:(m-1)){     x1[i+1]=(a*x1[i]+c)%%m}               m=1216     x=c()     int=c()     U=x1/m     for(i in 1:m){     x[i]=sqrt(9*U[i])-2     int[i]=exp(x[i]^2)*(x[i]+1)/(4+2*x[i])     }     I=9/m*sum(int)     I


title: Result from Step by step PCA implementation differs from `matplotlib.mlab.PCA()`, would be nice if someone can help me finding the source

text: I was reading this nice article and tried to implement this step by step guide in Python, and then I compared the results using the Python function from the matplotlib library: `matplotlib.mlab.PCA()`. However, I noticed that the results differ. Now, I am a little bit concerned if I have made a mistake in my implementation somewhere, or if the `PCA()` function takes a slightly different approach. would be cool to get some feedback on this very rough draft:   view iPython Notebook Thanks!


title: What is a relevant level for percentage deviance explained in a GLM?

text: I am trying to sort out the unique effect of various environmental predictors on species occurrence (presence/absence data). I have been running `glm` models in R with `family=binomial`. Most of my variables have very highly significant P values, but I have a large data set (~8000 data points), so this is maybe not so surprising. I have also been calculating the percentage deviance explained using the `BiodiversityR` package. Some of my very significant variables based on the P value have very low explained deviance (1% or so). Is there a rule of thumb for saying a predictor is not significant based on the deviance explained? Any other tests I should be including? I wanted to look at each variable on its own, and then test out specific interactions, but only if the variable was "important" in a model where it was the only predictor.


title: Estimating the parameters of a model, which method should I use?

text: I am trying to estimate a system of macroeconomic (simultaneous) equations, and I've learned about the 'existence' of various methods including Structural Equation Models, Simultaneous Equations Models, Generalised Methods of Moments, and Structural Vector Autoregression (not going into Bayesian methods). I know there is no short answer to this, but I would really appreciate if somebody could help me differentiate these methods. The objective of estimating the parameters is to understand the structure of the economy so that I may re-fit the model with a different set of data (for e.g. what happens to the output of the economy if the inflation target of the monetary authority is 4% instead of 0%). Thanks a lot!


title: Computing probabilities of letters in a word grid

text: I'm sure most people are familiar with word grid games like Boggle and the newer digital versions Scramble with Friends and Ruzzle. ![4x4 word grid](http://i.stack.imgur.com/Q2ePS.png) For anyone not familiar, the idea is to find words by using adjacent tiles. You start from any cell and try to spell a word by dragging up, down, left, right, or diagonal. The board doesn't wrap, and you can't reuse letters you've already selected. I'm trying to figure out the likelihood of a word appearing on a board given the likelihood of knowing how often individual letters appear. For example, if I know that the letter A appears 9.8% of the time, what is the probability of seeing the word AA? I know this is fairly simple, but it's been too long since college stats. (I have two hokey models, but I'd like to hear from the experts.) I could run a simulation of a million boards and come up with an empirical answer--in fact, someone has--but I'd rather understand why that is. In order to make things simpler, I'd like to ignore two rules that add to the complexity:   * We can ignore the constraints of how many times a letter can appear on a board. e.g., don't worry about whether there will be enough E's to make the word ELECTEE.   * We can ignore the fact that letters have to be adjacent. This means we don't have to figure out the likelihood that a letter is a neighbor of another letter that we need for the word So with that being said, and with the following probabilities   * A: 0.098   * E: 0.146   * R: 0.079   * S: 0.102   * T: 0.098   * C: 0.021 What is the likelihood of a random board containing the word SEE? What about SET? TEAR? What about the rarer CREATES? NOTE: I did search here already and I think this is a similar but different problem than this question: Probability of drawing a given word from a bag of letters in Scrabble


title: I have data that includes a record of each event and its time

text: I want to see if this other data I have predicts that this first event will occur. In other words, I have a data field with a record of farts, with columns for day, month, year, and then I have another data field with a record of meals, with columns for date month and year, and also a column for whether the eater is a bean eater. How do I analyze this? It would be easier if this data was in a easy time series format, with dates on the x axis and number of farts on the y axis, but the data only records when a fart occurred, not when a fart didn't occur. Thanks!


title: mvoutlier vs influence measures

text: I am exploring the mvoutlier package and comparing it against conventional influence measures such as Cook's Distance, Leverage, DFFITS etc. I have not been able to get my hands around comparative advantage of using one technique over the other. In your experience in which situations one is better than the other?


title: How do I report if linear regression and 2 way anova main interaction is providing me with different results for the same hypothesis?

text: I am doing a study on locus of control and coping ability moderated by social support. For the subscale of locus of control - internal belief, it was found that there were no significant prediction on coping ability when using the linear regression statistical analysis. However, when the internal belief subscale was categorized into high and low levels and tested with a 2 way anova on coping ability, it can be seen that there was a significant main effect. So, how should i report this contrasting find in my research study? It is confusing and I am stuck in the discussion section because of this contrast in results obtained. Please help and advice.


title: Variance covariance matrix for the inverse of a matrix

text: Suppose we have a matrix $\mathbf{A}=\begin{bmatrix}a_{11} & a_{12}\\\ a_{21} & a_{22} \end{bmatrix}$ and know its variance covariance $\left(4\times4\right)$ matrix. Then how the variance covariance matrix of $\mathbf{A}^{-1}$ can be obtained? I would highly appreciate your help. Thanks


title: Normality of residuals in a regression model with a categorical IV

text: I have ran a simple regression with a continuous response variable and a categorical explanatory variable (with 2 levels). I am currently checking that the model meets the assumptions of regression. I produced the following plot: ![enter image description here](http://i.stack.imgur.com/7s0f8.png) I'm aware that I need to check that the residuals are normally distributed. Do I need to check the distribution of residuals at each of the 2 levels of the explanatory variable? Or do I need to check the distribution of all residuals simultaneously?


title: Two sample t-test vs regression

text: I would like to examine the difference in forecast error between year 1 and 2. The descriptive statistics show that the mean of forecast error is higher in Year 2 but the median is lower. Under what circumstance would this happened? Also, when I test this in the following regression after controlling for other factors, $\beta_1$ is negative and significant, which means forecast error in Year 2 is significantly lower than Year 1. However, This result is inconsistent with the descriptive statistics result. How do I interpret that? Under what circumstance would the regression direction completely opposite to the t-test direction? $Error = \alpha + \beta_1Year + \text{Control variables}$ where: $Year$ is a dummy variable (year 2 is $1$, & year 1 is $0$).


title: Questions on standard deviation of a time series

text: If the standard deviation of a economic time series is approximately proportional to its level, that is, the standard deviation is well expressed as a percentage of the level of the series, then the standard deviation of the natural logarithm of the series is approximately constant. Could you tell me why? Thank you


title: Predictive models with large numbers of missing values in the features

text: I have been trying to train an algorithm to predict if an account will close or not using thousands of data points and many features. I am using data from the month before the account closed but the problem is that accounts have been around for different amounts of time. So, whereas for one account I might only have performance data up to 1 year, for another account I might have 1 month, 3 month, 1 year, 3 year, 5 year and even 10 year. There were 94 features but I cut it down to 19 to start playing around with it. I am looking for some help on how I can build an algorithm that incorporates the different amounts of data for each account. I was thinking about using a neural network because I have always been interested in them but I am open to any suggestions. Basically, then, I have many missing values in my features. If I simply omit observations with many missing values, my data set becomes far too small to be useful. Is there a standard way of handling this type of missing data problem, or a particular algorithm or model than handles it well? ![enter image description here](http://i.stack.imgur.com/y1VVj.jpg)


title: Interpreting Augmented Dickey-Fuller in R

text: I know that this has been discussed before, but those discussions did not really answer my questions. I know how the ADF test works, but I am having trouble interpreting the output for the three options using the `ur.df` function in R (package: `urca`). Could someone walk me through the interpretations? More specifically, what are `tau1`, `tau2`, `phi1`, `phi2`, and `phi3`?               summary(ur.df(tcm.ts, type="none",selectlags="BIC"))     summary(ur.df(tcm.ts, type="drift",selectlags="BIC"))     summary(ur.df(tcm.ts, type="trend",selectlags="BIC"))          > summary(ur.df(tcm.ts, type="none",selectlags="BIC"))          ###############################################      # Augmented Dickey-Fuller Test Unit Root Test #      ###############################################           Test regression none                Call:     lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)          Residuals:         Min      1Q  Median      3Q     Max      -95.199 -23.380  -6.608  26.885  86.560           Coefficients:                Estimate Std. Error t value Pr(>|t|)        z.lag.1     0.04398    0.01205   3.650  0.00183 **     z.diff.lag -0.03722    0.24417  -0.152  0.88053        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 45 on 18 degrees of freedom     Multiple R-squared:  0.7091,    Adjusted R-squared:  0.6768      F-statistic: 21.94 on 2 and 18 DF,  p-value: 1.492e-05               Value of test-statistic is: 3.6495           Critical values for test statistics:            1pct  5pct 10pct     tau1 -2.66 -1.95  -1.6          > summary(ur.df(tcm.ts, type="drift",selectlags="BIC"))          ###############################################      # Augmented Dickey-Fuller Test Unit Root Test #      ###############################################           Test regression drift                Call:     lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)          Residuals:         Min      1Q  Median      3Q     Max      -69.366 -24.625  -3.018  34.165  82.227           Coefficients:                  Estimate Std. Error t value Pr(>|t|)       (Intercept) -75.21181   45.89715  -1.639   0.1196       z.lag.1       0.09756    0.03467   2.814   0.0119 *     z.diff.lag   -0.20396    0.25469  -0.801   0.4343       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 43.03 on 17 degrees of freedom     Multiple R-squared:  0.3596,    Adjusted R-squared:  0.2843      F-statistic: 4.773 on 2 and 17 DF,  p-value: 0.02264               Value of test-statistic is: 2.814 8.6258           Critical values for test statistics:            1pct  5pct 10pct     tau2 -3.75 -3.00 -2.63     phi1  7.88  5.18  4.12          > summary(ur.df(tcm.ts, type="trend",selectlags="BIC"))          ###############################################      # Augmented Dickey-Fuller Test Unit Root Test #      ###############################################           Test regression trend                Call:     lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)          Residuals:        Min     1Q Median     3Q    Max      -86.46 -14.84   5.56  20.87  70.29           Coefficients:                 Estimate Std. Error t value Pr(>|t|)       (Intercept)  91.4039    92.1407   0.992   0.3360       z.lag.1      -0.1127     0.1082  -1.042   0.3129       tt           13.0810     6.4318   2.034   0.0589 .     z.diff.lag   -0.1287     0.2369  -0.543   0.5946       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 39.54 on 16 degrees of freedom     Multiple R-squared:  0.4911,    Adjusted R-squared:  0.3957      F-statistic: 5.148 on 3 and 16 DF,  p-value: 0.01109               Value of test-statistic is: -1.042 8.1902 6.7579           Critical values for test statistics:            1pct  5pct 10pct     tau3 -4.38 -3.60 -3.24     phi2  8.21  5.68  4.67     phi3 10.61  7.24  5.91


title: Bodyfat, abdomen and biceps: How to interpret this data?

text: > > With the help of the abdominal girth and the girth of the biceps we want > to prophesy the bodyfat. To do so we collected the data of 252 men (the > first 40 data are shown below). First do the linear regression of the model > bodyfat/ abdomen, then the liear regression of the model bodyfat/ biceps and > then the multiple regression of the model bodyfat/ (abdomen + biceps). > Compare and interpret the results! Below is a part of the data in order to give you an impression what the data look like. With R I did the two simple linear regressions and the one multiple regression and drew it (see below). **Model 1** _(bodyfat/ abdomen)_ For the linear model I got the two estimated coefficients $$ \hat{\theta}_1\approx-35.197,~~~\hat{\theta}_2\approx0.585. $$ With this I drew the line into the data that you see in the first picture. **Model 2** _(bodyfat/ biceps)_ Here, I got the two estimated coefficients $$ \hat{\theta}_1\approx-21.882,~~~\hat{\theta}_2\approx1.265. $$ Again I drew the line (s. second picture below). **Model 3** _(bodyfat/ (abdomen + biceps))_ Here, three coefficients had to be determined, I got $$ \hat{\theta}_1\approx -30.684,~~~\hat{\theta}_2\approx 0.645,~~~\hat{\theta}_3\approx -0.311. $$ Now the task is to compare and to interpret. But to be honest, I have some difficulties to read something in this... or to interpret something. Maybe you can help me to compare and to interpret? What can be seen from all this? * * * ![data](http://i.stack.imgur.com/SQyUv.png) ![graphics and regressions](http://i.stack.imgur.com/9SwLZ.png) **EDIT** The $R^2$-value of the first model is $R^2=0.6621$. The $R^2$-value of the second model is $R^2=0.2431$. The $R^2$-value of the third model is $R^2=0.6699$. So between the abdomen and the bodyfat there is a rather good linear connection, but the linear connection between the biceps and the bodyfat seems to be rather small. Put abdomen and biceos together, there is a better linear connection which seems to be clear, because the bodyfat depends on very much factors (such as food, moving, gene, ...) and the more factors one puts into consideration, the better is the connection.


title: how to fit a zero inflated NB distribution with large mu and sigma by fitdist + gamlss functions

text: I've been trying to figure out how to use gamlss.famliy Zero inflated functions (e.g. ZIP) on the data. A descriptive stat on my data is as follows;               Min. 1st Qu.  Median    Mean   3rd Qu.    Max.        0     3       45      13520    630    469000       A result of `descdist` ;               min:  0   max:  469000      median:  45      mean:  13522.79      estimated sd:  66346.55      estimated skewness:  5.996958      estimated kurtosis:  39.37641       I could run `fitdist(W14,"nbinom")` as follows; ![enter image description here](http://i.stack.imgur.com/CkUwi.png) However, when I run `fitdist(W14, "ZIP2", start=list(mu=x, sigma=x))`, and I get a usual error message;               sigma must be between 0 and 1      Error in fitdist(W14, "ZIP2", start = list(mu = xx, sigma = xx)) :      the function mle failed to estimate the parameters,      with the error code 100      I've tried the different values of mu and sigma but none has worked so far. I really do appreciate any help on this.


title: Comparing proportions in two samples

text: (x-posted to Statalist) I have two (probably) independent samples taken from the same population under two different surveys conducted a few years apart. I'd like to check whether the frequencies for a particular indicator match between the samples. I don't expect these proportions to change much over time (though I can't rule it out). The indicator uses multiple categories which are mutually exclusive and not ordinal, e.g.: What is your car's color?   1. Red    2. Green   3. Black   4. Blue etc. I want to check that survey A captured, basically, the same proportions of red/green/blue that survey B did. I _think_ I should be using a Chi-square goodness of fit for this (`csgof`, in Stata); except that I'm comparing two samples, rather than the sample versus the population. Anyone know what other statistical test I should be using? I'm a bit stuck, and am considering just eyeballing it: using `svy: tab` (in Stata) on both and just comparing the confidence intervals. But I feel like this is a crude way of doing things. Thanks very much! a


title: Clustering communication patterns to detect multiple identities

text: I have a data set of communication patterns between chatting agents. Each agent can have multiple profiles or identities. I am interested in developing a way to investigate the similarity between various profiles. The problem is that there is no real indicator whether two profiles belong to the same agent. The set of features describes the frequency of messages through the hours of the day, days of the week, months, and years. The aim is to learn the significant features that indicate similar profiles that belong to one agent. The proposed approach: considering each data item to be the communication pattern of one agent over a year. This way, the data set would have many data items that represent the communication patterns of the same agent across different years. We apply clustering with a number of clusters equal to the number of agents--not profiles. Ideally, if the clustering algorithm clusters the data points that belong to one agent in one cluster then we have a features model that we could use to detect profiles of the same agent. Of course this proposal assumes that the agent would have the same pattern in the different profiles, which is not expected to be the case. But among the data set of the all agents, this way may give a correct prediction with a certain level of accuracy. Does this approach sound reasonable? Is there a better way to tackle this issue? Thank you!


title: R web interface suggestion

text: hi i was looking for a sophisticated web based interface of R which can be used for couple of analytics application? any suggestion please..rweb is quite simple


title: Determining the significance of a set ordering

text: First some preamble: Let A be a set of unordered elements, and let A' and A'' be subsets of A, where |A'| = |A''| = k. It’s given that as k approaches |A|, for every x in A' there is an increased probability of finding x in A''. That is, as k approaches |A|, |A' ∩ A''| approaches |A|. On to the issue at hand: Let’s say we have two different graphs G1 and G2, each with n vertices labeled 1 through n, and each with a set of non-random edges.* Let B be the set of vertices from G1, and C be the set of vertices from G2. Order B and C by vertex degree. Let B' and C' be subsets of B and C respectively, where |B'| = |C'| = k, and where subset elements are taken from their respective supersets from left to right (see example below), such that as k approaches n, B' and C' will tend to include vertices of lesser and lesser degree. Here’s a general example of the subset selection process: B = {1,4,2,5,3} C = {4,2,3,1,5} k = 1 -> B' = {1}, C' = {4} k = 2 -> B' = {1,4}, C' = {4,2} k = 3 -> B' = {1,4,2}, C' = {4,2,3} k = 4 -> B' = {1,4,2,5}, C' = {4,2,3,1} k = n -> B' = {1,4,2,5,3}, C' = {4,2,3,1,5} If I were to plot k vs |B' ∩ C'|, the shape of the curve would be attributable both to the effect of increasing k, as described in my preamble, but ALSO, and especially, to the effect of the ordering. In this case, the ordering by degree. How do I determine the significance of this ordering? *By non-random edges I mean that the connectivity of the respective graphs is determined by some other factor. In my real-world case, these graphs represent brain connectivity at different intervals of time. Connected nodes represent highly correlated activity between neurons.


title: Is it valid to use Hartigans' dip test to reject uni-modal null hypothesis with large N?

text: first time questioner, so please be gentle :) I have two distributions of data from a simulation. By eye, one looks like it might be bimodal, one not. I copy them below ![enter image description here](http://i.stack.imgur.com/Kl6fa.png) **A: parameter value 1 :** by eye possibly bi / multi modal ![enter image description here](http://i.stack.imgur.com/nF8Fc.png) **B: parameter value 2 :** by eye unimodal I thought I would use Hartigans' dip test to test H1:not uni-modal vs H0:uni- modal. The problem is that the test statistic suggests that I reject the null hypothesis of unimodality for _both_ distributions with values well under the suggested 0.05 threshold. In fact, the test statistic for the distribution that looks more unimodal is lower than that for the one that looks potentially multi-modal (dist A: D = 0.00814; dist B: D = 0.00340) I think what I'm seeing is the effect of fairly large N (N=5000), so the sample size lends statistical power to the test. But inspection of the histograms suggests this is not valid. Is there some way to formalise discussion of whether it is valid to reject the null hypothesis based on this test? I've read a few posts on here (Test for bimodal distribution and @whuber 's suggested search ). I also found this elsewhere, but the answer is a bit generic - basically saying that with large N you're likely to find lots of tests significant, which I already suspect is the case here. I realise that some consideration of causal mechanisms for uni / non-uni modal outcomes could aid the discussion, but want to understand the statistical tests, too. I'd like some advice on 1) Am I interpreting the test statistic correctly (i.e. D < 0.05 evidence to reject uni-modality)? 2) Is there a way to determine whether large N is giving the test undue statistical power?


title: Do we need an unbiased estimator of the variance?

text: "Although it is nice to have an unbiased estimator of the variance, we do not really need it to understand the relation between our independent variable and our dependent variable. Why?" I think I kind of get it, but he wants a specific answer, and he picks on me...


title: Overfitting a linear Linear Discriminant Function

text: I am estimating a Linear Discriminant function with 250 input variables over 4000 data records. Should I consider feature selection, am I over fitting the model? How do I know when feature selection should be considered? I'm using a Linear Programming Discriminant Analysis model proposed by Stam and Ragsdale