title: Is there a way to perform SVD in a sequential manner?

text: My neurology experiment has a spike detector outputting `40` sample long spike waveforms. I'm using a dictionary method for sorting the spikes in real time. To obtain the dictionary at start of program, I want to use all the spikes detected in the first `30` seconds and perform an `SVD` of the `Nx40` matrix, where `N` is the number of spikes detected in the first `30` seconds. Is there an algorithm which can perform this `Nx40` `SVD` in a sequential or iterative method or perhaps in smaller batches. The expected spiking rate means my program cannot store all the 30s worth of spike waveforms (rows of the matrix) in memory.


title: Understanding vec2var conversion in R

text: I'm using Bernhard Pfaff's packages {urca} and {vars} to analyze 3 time series. Each is I(1) and cointegrated with $r =2$ cointegrating relationships. The _vec2var_ () command should make the conversion from the **ca.jo** object (my VECM) to its VAR representation in levels. But I do not understand the output from _vec2var_ (), nor how to perform VAR analyses on this output (i.e., send the data to the _VAR_ () function to get a _varest_ object). That is, after converting from VECM to VAR what do I send to the _VAR_ () function (which performs the standard VAR estimates)? Should it be:               var.form$A%*%var.form$datamat[ , c(6:11)] + var.form$deterministic%*%var.form$datamat[ , c(4:5)]      where dimensions in the code above are because I have $K=3$ variables, $p=2$ lags, a constant term and a deterministic trend. But the formulation from _vec2var_ () [on which I've based the above R code] is something like \begin{equation} Y_t = A_1 Y_{t-1} + A_2 Y_{t-2} + \Phi D_t + \mu + \epsilon_t \end{equation} That is, a basic VAR(2) with a deterministic time trend and constant. But this is what I'm confused about: The VAR form of the long-run VECM (equation 4.8a, in Pfaff's _Analysis of Integrated and Co-Integrated Time Series with R. useR!_ , and elsewhere _e.g._ Lütkepohl eq'n 7.1.1) is: \begin{equation} \Delta y_t = \Gamma_1 \cdot \Delta y_{t-1} + \Gamma_2 \cdot \Delta y_{t-2} + \Pi \cdot y_{t-2} + \Phi \cdot D_t + \mu + \epsilon_t \end{equation} To get the proper VAR level representation do I have to take the data from _vec2var_ () and manipulate it myself? (i.e. find the difference series and then multiply by the appropriate cointegration matrix)? If this is the case then **which** [in R code please!] coefficient matrices should multiply which vectors? I feel like there's a disconnect from the theory and mathematical formalism and these R packages... but more likely I'm missing something pretty basic. Here's the basic R code set up:               library(urca)      library(vars)           vecm <- ca.jo(mydata, ecdet = "trend", type = "trace",  spec = "longrun")      var.form <- vec2var(vecm, r= 2)      # The output data I get can be seen from names(var.form), but the key parts are:           var.form$y          # my original data      var.form$datamat    # original data organized with the three variables (starting at t = 2), the constant term, trend term and then six columns with two lagged levels at t=1 and t=0.           # the coefficients for the regressors on the endogenous variables (A) and deterministic components are:          var.form$A             # lagged variable coefficients (levels only, I think)           var.form$deterministic     # constant and time trend coefficients


title: How to compute significant interaction estimates when main effect is not significant?

text: I have a linear model of a dependent variable, $y$, with two predictor variables, year and site, and their interaction, with year being numeric and site categorical. The main effect of year is not significantly different from zero, with an estimated value of 0.02312 for the $y$ on year slope. Some of the year by site interactions are significant. The summary of the linear model in R gives estimates for the interaction terms as deviations from the main effect of year. If, for example, the year by site 1 estimate is reported as .416 and is significant, then to compute the $y$ on year slope within site 1, should I compute 0 + .416 or 0.02312 + .416?


title: How to deal with repeated measurements in the same condition of a factorial experiment?

text: I am in Psychology and trying to explore the utility of mixed modeling for analyzing my repeated-measures data in a factorial experiment. The primary reason for using mixed models is that I would like to avoid the common practice of averaging data collected in the same experimental condition. My understanding is that it's typically required for repeated-measures ANOVA that there is only one observation per condition per subject. What if you have several replications of the same condition for the same subject? To be more concrete, I have two conditions a between-subjects factor A (2 levels) and a within-subjects factor B (3 levels). There are 4 repetitions of each level of B for a total of 12 randomized order trials per subject. Usually, I would simply average across these 4 to get an estimate for the performance of the subject in the condition and then run, but it seems that this way I'm throwing valuable information about variability. How to deal with such data using mixed modeling in R (I've been using lmer function). Maybe including trial number as another variable would work? I tried including trial # as a random factor together with the subject, but its estimated variance is very low compared to error and subject.


title: Joint likelihood of two models

text: I am estimating a hurdle model, which consists of two data-generating processes: one for the zero/positive relationship, and a different one for the positive numbers above zero. A binary logit model is often used for the first part, and a zero-truncated poisson or negative binomial model for the second. In R, the `hurdle()` command in the pscl package jointly estimates the models together, and returns a single AIC and log-likelihood value. I need to incorporate some random effects, and therefore am using the `glmmADMB()` command in the package of the same name. This package does not have a single hurdle model command, so I must estimate them separately (first binary, then poisson). My question: how should one calculate the joint log-likelihood of two separately estimated models? I assume this is possible, given that `hurdle` seems to do just this.


title: Time series forecasting lookback windows -- sliding or growing?

text: Are there any good reasons to prefer a sliding model training window to a growing window in online time series forecasting (or vice versa)? I'm particularly referring to financial time series. I would intuitively think a sliding window should perform worse -- out of sample-- as it is has more potential for over-fitting specific sample window characteristics, but some of the empirical results I've seen are counter to this. Also, given that a sliding window is preferred by some, what would your approach be to determine the look-back length (any good reasons to prefer one over another, aside from pure heuristics)? Although I didn't specify a model, an example might be ARIMA. EDIT: I should add there there is a related blog post by By. Rob Hyndman, with what he dubbed 'time-series' cross validation. While it does cover the concepts described, it doesn't give much of a formal reason about why one method might be preferential over the other, nor any ideas about an optimal look-back window parameter.


title: Smoothing a time series of ratios

text: I have a time series of proportions, $x_t = \frac{a_t}{b_t}$ i.e. $x_1 = 2 / 30, x_2 = 1/10$, ... I want to smooth $x_t$. Should I apply a smoothing function directly to $x_t$, or should I smooth $a_t$ and $b_t$ individually, then recompute the proportions? My gut tells me you should smooth the individual time series, and then recompute. Suppose each of the three time series have noise (this is why we would be interested in smoothing to begin with), then maybe $a_t = f(t) + \epsilon_a$, $b_t = g(t) + \epsilon_b$. Then the proportions $x_t = \frac{f(t) + \epsilon_a}{g(t) + \epsilon_b}$, which could be smoothed by taking $\hat{x_t} = \frac{\hat{f(t)}}{\hat{g(t)}}$


title: Statistics usage for customer event

text: We are preparing one event for our customers and I need to count numbers we need to prepare for this event (I think my question belongs to statistics and it isn't a homework). Inputs:   * We have 26 ASCII letters (ASCII codes 65-90) from which customers choose. When somebody choose a letter, he doesn't give it back.   * We have 500 playing cards (= 500 customers can compete) with three empty boxes for letters.   * In one draw customer gets one letter per one draw lot   * The winner is somebody who has three particular letters, e.g. "A","L","Q".  Question: How many letters we need to prepare if we have prices for 100 winners?


title: What are i.i.d. random variables?

text: How would you go about explaining i.i.d to non techncial people?


title: Using the Bayes Theorem?

text: > A certain town has two taxi companies, the Green Taxi Co (cars coloured > green) and the Blue Taxi Co (cars coloured blue). 10% of taxis are the Green > and 90% are the Blue. There was an accident on a dark evening, and the > witness claimed a Green cab was involved. On future testing, they discovered > that under the conditions of that evening, there was an 80% chance of > correct identification of the cab colour (regardless of colour) and 20% > chance of erring in colour identification. (assume neither company is more > accident prone than the other) > > a. Before actually hearing the testimony of this witness, what is the > probability that the witness will claim that it was a Blue Cab involved in > the accident that night? > > b. If we now have the testimony of the witness, who declares that the cab > was Green, what is the probability that it actually was Green? I used a probability tree to solve this question but apparently, it was actually related to Bayes theorem. How I did (a) was use the probability tree and got $P(Witness~claims~Blue) = (0.9)(0.8)+(0.1)(0.2)=0.74$. For (b), $P(Actually~green) = (0.1)(0.8)+(0.1)(0.2)=0.10$. However, the answer for (b) was $\frac{.08}{.26}$ and I was supposed to use the Bayes Theorem. Can someone please explain why the Bayes is to be used instead of a sum of probabilities?


title: Is a larger beta weight a better predictor than a high t-statistic?

text: Say there are two variables. Variable A has a higher $\beta$ weight than variable B but variable B has a higher t-statistic than variable A. Can either of these measurements, standardized $\beta$ or t-value, offer better predictive power over the other? Must t-value and $\beta$ be taken into consideration together or not?


title: Can I do a correlation between a likert scale and a dichotomous scale with the same list of items

text: I am trying to do a correlation between a 4 level likert scale with 11 items and a dichotomous present-not present scale of the same 11 items. Can I do this?


title: Classification of states in Markov Chain

text: **Question** Consider the following transition matrix:                P=       0     0     1     0     0     0       0     0     0     0     0     1      0     0     0     0     1     0     1/4   1/4    0    1/2    0     0      1     0     0     0     0     0      0    1/3    0     0     0    2/3      **a)** Which states are transient? **b)** Which states are recurrent? **c)** Identify all closed sets of states. **d)** Is this chain ergodic? Dear friends , I have thought the states are as follows, {1,5} {0,2,4} recurrent since they communicate with each other {3} transient and since state 3 does not communicate with other states, it is not an ergodic Mc Am I correct? Thank you..


title: Explain what a kernel is, in plain English

text: there are several distinct usages: kernel density estimation   kernel trick   kernel smoothing Please explain what the "kernel" in them means, in plain English, in your own words


title: Exact probability and asymptotic approximation

text: From Wikipedia > A large amount of computing resources is required to compute **exact > probabilities** for the Kruskal-Wallis test. Existing software only provides > **exact probabilities** for sample sizes less than about 30 participants. > These software programs rely on **asymptotic approximation** for larger > sample sizes. **Exact probability values** for larger sample sizes are > available. Spurrier (2003) published exact probability tables for samples as > large as 45 participants. Meyer and Seaman (2006) produced exact probability > distributions for samples as large as 105 participants. Critical value > tables and exact probabilities from Meyer and Seaman are available for > download at http://faculty.virginia.edu/kruskal-wallis/. A paper describing > their work may also be found there. I think exact probabilities are the accurate probabilities, but why does it say that "Existing software only provides **exact probabilities** for sample sizes less than about 30 participants. These software programs rely on **asymptotic approximation** for larger sample sizes"? Is it wrong? Thanks and regards!


title: How do you compare two Gaussian Processes?

text: Kullback-Leibler divergence is a metric to compare two probability density functions, but what metric is used to compare two GP's $X$ and $Y$?


title: Mixed model multiple comparisons for interaction between continuous and categorical predictor

text: I would like to use `lme4` to fit a mixed effects regression and `multcomp` to compute the pairwise comparisons. I have a complex data set with multiple continuous and categorical predictors, but my question can be demonstrated using the built-in `ChickWeight` data set as an example:               m <- lmer(weight ~ Time * Diet + (1 | Chick), data=ChickWeight, REML=F)      `Time` is continuous and `Diet` is categorical (4 levels) and there are multiple Chicks per diet. All the chicks started out at about the same weight, but their diets (may) affect their growth rate, so the `Diet` intercepts should be (more or less) the same, but the slopes might be different. I can get the pairwise comparisons for the intercept effect of `Diet` like this:               summary(glht(m, linfct=mcp(Diet = "Tukey")))      and, indeed, they are not significantly different, but how can I do the analogous test for the `Time:Diet` effect? Just putting the interaction term into `mcp` produces an error:               summary(glht(m, linfct=mcp('Time:Diet' = "Tukey")))     Error in summary(glht(m, linfct = mcp(`Time:Diet` = "Tukey"))) :        error in evaluating the argument 'object' in selecting a method for function      'summary': Error in mcp2matrix(model, linfct = linfct) :      Variable(s) ‘Time:Diet’ have been specified in ‘linfct’ but cannot be found in ‘model’!


title: time series --- seasonal adjustment

text: I'm concerned to seasonal adjustment procedure and want to know the criteria for this purpose can anyone please give me the answer of the following question. what should be the criteria for seasonal adjustment whether to do or not to do?


title: Examine neg. correlated relationships

text: I have some environmental data I am examining (samples from water of a lake). After next generation sequencing and some computation I get some OTUs (Operational Taxonomic Units) saying which bacteria are in which environmental samples (150 samples). First I filtered for the OTUs which have a high abundance in general and occur at least in five samples. Then I calculated spearman correlation-coefficients between all the OTUs. With the positive coefficients I did then a Markov clustering which gave me some interesting results. The question is now how to use the negative correlations (In combination with the pos. correlations or alone). I want to investigate in the question, which bacteria have a strong "avoidance" (from competition, exclusion by environment, predation, or whatever). What method would you suggest? Classical clustering methods always use the idea of similarity or distance, but this does not make sense here. For every sample I have additional information like temperature, depth, date, oxygen and many more.


title: Generating samples of correlated normally distributed variables with some of the variables pre-selected

text: I would greatly appreciate any of you who could help me with this challenge. I am going to state the problem in sequential order, so as to make it clear:   1. I have $n$ normally distributed random variables.   2. I have statistical data for each random variable (thus, their mean and variance).   3. I have historical data for each of the random variables, and from these data I can calculate the correlation between the random variables. Now, assuming I know the value of one of the variables, how do I select the values of the other variables, such that all the $n$ values in the resulting vector are correlated accordingly. If that is possible, how do I extend this to a situation where I have two or more known values? For example, I know two or more values out of the $n$ variables (they satisfy the correlation constraint), and want to pick the others so that they also satisfy the correlation and their respective distribution functions. I am not well-versed in statistics, so I would appreciate explanations that are as simple to understand as possible.


title: How to extract "average polyline" from a set of polylines

text: EDIT: As pointed, there is already a similar question "buried" inside a larger-scope one. I'll reproduce the relevant part here: > If I have multiple recordings of the same routes are there any valid > techniques to combine them to get closer to the real route? The relevant part of the chosen answer, in the context of my question, is this: > With multiple recordings, you can create a 2D kernel smooth of each, sum the > smooths, and subject that to a topographic analysis to look for a > "ridgeline." You won't get a single connected line, usually, but you can > often patch the ridges together into a continuous path. People have used > this method to average hurricane tracks, for example. * * * I have a set of polylines in bidimensional space (X,Y) representing object trajectories. Every polyline in the set is a measured track along the same prescribed path, most of them (but not all) starting around the same point, and ending around the same other point. I would like to statistically calculate the idealized path ("average polyline") from the set of polylines, in the same way I would manually manually trace it in a drawing program, except I would like this to happen without manual intervention. A sample image of a part of the polyline set (cropped) is below. It seems to me that tere is an "obvious" place where I would draw the average line with the mouse, but I can't figure which algorithm I could use. Important notes:   1. I know I could use the nodes for calculations, but some segments might be arbitrarily long, thus not having nearby nodes that would help in the calculations. That is, I would like to consider the line segments themselves for the calculation, not just the nodes.   2. There might be parts of a polyline behaving like outliers (see image). I think that, given the majority of polylines is well behaved, this outlier part would easily be "absorbed" by the well-behaved ones nearby.   3. Although the original data is time-stamped, I believe no temporal relation should be necessary to infer the spatial result I am after, isn't it?   4. Yes, this is GPS data (although this is incidental to the problem formulation, which has been purposefully abstracted).  ![enter image description here](http://i.stack.imgur.com/FoIjE.png)


title: Does GLMM or LMM already give us planned contrasts?

text: I use GLMM and LMM to find significant factors for my data. And could I ask if the results from these two models already give me planned contrasts? Or do I need to run them separately?


title: Fit Gaussian random field for spatial data

text: I'm dealing with spatial data where the response variable is the gas concentration. In addition, I've the x,y-coordinate values, and another covariates. I'm thinking to fit a Gaussian random field model as follows: $$ Z(s) = X(s) \beta + \epsilon(s) $$ where $ [Z(s) : s \in R^2]$, $X(s)$ is a set of p covariates associated with each site s, $\beta$ is a p-dimensional vector of coefficients. I already estimated the covariance function for the error term, but now I'm having trouble on how to proceed. I'm using R, but so far, unfortunately, I've not come across any package where I can include covariates in the model. I'm still learning spatial statistics, so if you've any hint, I would truly appreciate it.


title: glmnet, categorical variable, group lasso?

text: I am using glmnet for LASSO. My data set contains several continuous variables and one categorical variable (it has four levels). I wondered if I could treat three dummy variables as other continuous variables. Should I use a type of group LASSO approach for the three dummies? I'm sorry if this is a dumb question.


title: What exactly is the actual coverage probability?

text: I am confused about the difference between the **nominal** and the **actual** coverage probability. Say we are trying to estimate the fraction of the current population of the US that has been diagnosed with cancer. We get a sample (e.g. 10,000 individuals) and then estimate this fraction ($f$) and a confidence interval around it $[a,b]$ with $\alpha$ (confidence level) $=0.95$. At this point, if I understood things correctly, the nominal coverage probability is $\alpha$, i.e. $0.95$. But what about the **actual** coverage probability? To me, there is only one true value of $f$, e.g. say $5\%$, and it doesn't make sense to talk about the _actual_ probability that $f$ is within the interval $[a,b]$. It either is or it isn't within $[a,b]$. What's wrong with this reasoning? And how would one obtain the actual coverage probability in practice?


title: What's the practical difference between the Johansen vs Engle-Granger tests for cointegration?

text: For the two-variable case, what are the practical differences between using the Engle-Granger procedure versus the Johansen test for cointegration? Is one universally more powerful than the other? Will one give more false positives or false negatives than the other? Should Johansen always be preferred?


title: Reflect & Transform Negatively Skewed Data

text: I have data that are non-normal and (strongly) negative skewed. The data also have high kurtosis and outliers. There appears to be a variety of options for transformation, but I cannot find a source that helps me determine the best option (or why to choose one over the other). In this case, should I   1. reflect and transform    2. perform an exponential transformation    3. look into Box-Cox?  And how do I know which one is the right choice for my data?


title: How to fit a logistic regression for 1 dependent variable and 1 qualitative variable measured twice

text: I am struggling to fit a simple logistic regression for one dependent value (group) by one independent qualitative variable (dilat) measured twice independently (rater). I try many solutions and think according http://www.ats.ucla.edu/stat/mult_pkg/whatstat/ that the solution is a Mixed Effects Logistic Regression.               glmer_dilat<-glmer(group ~ dilat + (1 | rater), data = ex, family = binomial)     summary(glmer_dilat)                    Generalized linear mixed model fit by the Laplace approximation      Formula: group ~ dilat + (1 | rater)         Data: ex         AIC   BIC logLik deviance      105.5 112.5 -49.74    99.48     Random effects:      Groups Name        Variance Std.Dev.      rater  (Intercept)  0        0           Number of obs: 76, groups: rater, 2          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)       (Intercept)   0.8473     0.4880   1.736   0.0825 .     dilat        -1.2827     0.5594  -2.293   0.0219 *      But the result is the same without !               summary(glm(group ~ dilat, data =ex, family = binomial))          glm(formula = group ~ dilat, family = binomial, data = ex)          Deviance Residuals:         Min      1Q  Median      3Q     Max       -1.552  -0.999  -0.999   1.367   1.367            Coefficients:                 Estimate Std. Error z value Pr(>|z|)       (Intercept)   0.8473     0.4879   1.736   0.0825 .     dilat        -1.2826     0.5594  -2.293   0.0219 *      What is the solution? please find my data set here after applying a dput command to it.               structure(list(id = structure(c(38L, 11L, 6L, 28L, 3L, 30L, 39L, 4L, 8L, 12L, 32L, 29L, 34L, 35L, 33L, 16L, 27L, 5L, 36L, 10L, 9L, 14L, 1L, 13L, 31L, 2L, 17L, 7L, 19L, 20L, 18L, 21L, 23L, 15L, 24L, 25L, 37L, 26L, 38L, 11L, 6L, 28L, 3L, 30L, 39L, 4L, 8L, 12L, 32L, 29L, 34L, 35L, 33L, 16L, 27L, 5L, 36L, 10L, 9L, 14L, 1L, 13L, 31L, 2L, 17L, 7L, 19L, 20L, 18L, 21L, 22L, 23L, 15L, 24L, 37L, 26L), .Label = c("1038835", "2025267", "2053954", "3031612", "40004760", "40014515", "40040532", "40092413", "40101857", "40105328", "4016213", "40187296", "40203950", "40260642", "40269263", "40300349", "40308059", "40327146", "40333651", "40364468", "40435267", "40440293", "40443920", "40485124", "40609779", "40628741", "40662695", "5025220", "E9701737", "M/377313", "qsc22913", "QSC29371", "QSC43884", "QSC62220", "QSC75555", "QSC92652", "QSD01289", "QSD02237", "U/FY0296" ), class = "factor"), group = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), rater = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), dilat = c(1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L), midbrain_atroph = c(1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), quadrigemi_atroph = c(1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), hum_sig = c(1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), flower_sig = c(0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), fp_atroph = c(0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), scp_atroph = c(0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("id", "group", "rater", "dilat", "midbrain_atroph", "quadrigemi_atroph", "hum_sig", "flower_sig", "fp_atroph", "scp_atroph"), class = "data.frame", row.names = c(NA, -76L))


title: Should I use raw scores or profile scores for variables in a regression model?

text: I'm trying to build a regression model to predict subjective memory reports based on the following predictors: objective memory test, illness perceptions, verbal fluency, and mood. I also want to control for age, years of education and IQ by putting them into the regression model too. However, there's a problem. The objective memory memory test (The Rivermead Behavioural Memory Test) gives you both raw scores and profile scores for the various subscales. Some, but not all, of the profile scores are converted based on IQ and age. Also, the IQ estimate (The Test of Premorbid Functioning) converts raw scores to an IQ estimate, either just with the raw score, or with the raw score AND demographic variables of age, gender and years of education. The IQ estimate that includes demographic information is far more accurate than without. If I use profile scores (with demographics) then aren't I controlling for them twice? Should I use raw scores or profile scores? And what, and how should I enter these into the regression model? Sorry if the above doesn't make sense. I'm not sure I understand it myself. I'm really struggling to get my head around this! Thanks!


title: Comparing p-values from multiple permutation (Mantel) tests on variations of the same data set

text: I am comparing the relationship between genetic and geographic distance of individuals in a wild animal population. The hypothesis is that individuals with higher genetic relatedness establish home ranges closer to one another (shorter geographic distance between them) than to individuals that are less related (ie: geographic distance and relatedness are inversely correlated). Because these animals have homeranges that overlap and were sampled multiple times over the course of the study, I initially calculated a center point for each individual's home range and used that as the geographic location of the individual. I then ran a Mantel test (and the modification of the Mantel test called the Tau Kr test) and compared my correlation coefficient to the results of 1000 permutations of the data. I obtained a significant result (P<0.05). However, I wanted to see how sensitive the analysis is to where you sample the individual. Therefore, I re-ran my Mantel tests but using every possible combination of sampling locations for each individual (rather than taking the center point). I now have a distribution of p-values and test statistics for all of these tests and in 14% of cases the tests were significant (P<0.05) for the Mantel Test and in 52% of cases using the Tau Kr Test. Is there a way to evaluate if the initial hypothesis is true with these results? Intuitively to me, I think that if there was no relationship between geographic and genetic distance, even 14% of tests should not be significant. But I am not really sure what to compare this distribution of p-values to - what would the null model look like? It seems to be that I am detecting a real pattern even though it is highly dependent on sampling location, but I cannot figure out however how to analyze this data set now to determine whether the signal is significant. Many thanks for your help!


title: What is discriminative partition matching?

text: In the following bioinformatics paper, "Quantifying environmental adaptation of metabolic pathways in metagenomics", Gianoulis et al. employ the use of two tools to detect multivariate relationships between environmental features and microbiomic features:   * Regularized Canonical Correlation Analysis   * Discriminative Partition Matching (DPM) The paper references several statistics papers and books, including a book by R. Wichern, but it unclear which of these is the reference for DPM. A Google search for "Discriminative Partition Matching" pulls up several applied papers, including the Gianoulis paper, but no direct expositions of the method. Where can I find information on "Discriminative Partition Matching"?


title: R: Post hoc tests in ANCOVA with factorial design

text: Does anyone know how to do post hoc tests in ANCOVA with a factorial design? I have two vectors consisting of 23 baseline values (covariate) and 23 values after treatment (independent variable) and I have two factors with both two levels. I created an ANCOVA model and calculated the adjusted means, standard errors and confidence intervals. Example:               library(effects)          baseline = c(0.7672,1.846,0.6487,0.4517,0.5599,0.2255,0.5946,1.435,0.5374,0.4901,1.258,0.5445,1.078,1.142,0.5,1.044,0.7824,1.059,0.6802,0.8003,0.5547,1.003,0.9213)     after_treatment = c(0.4222,1.442,0.8436,0.5544,0.8818,0.08789,0.6291,1.23,0.4093,0.7828,-0.04061,0.8686,0.8525,0.8036,0.3758,0.8531,0.2897,0.8127,1.213,0.05276,0.7364,1.001,0.8974)          age = factor(c(rep(c("Young","Old"),11),"Young"))      treatment = factor(c(rep("Drug",12),rep("Placebo",11)))          ANC = aov(after_treatment ~ baseline + treatment*age)          effect_treatage = effect("treatment*age",ANC)     data.frame(effect_treatage)            treatment   age       fit        se     lower     upper     1      Drug   Old 0.8232137 0.1455190 0.5174897 1.1289377     2   Placebo   Old 0.6168641 0.1643178 0.2716452 0.9620831     3      Drug Young 0.5689036 0.1469175 0.2602413 0.8775659     4   Placebo Young 0.7603360 0.1462715 0.4530309 1.0676410      Now I want test if there is a difference between the adjusted means of Young-Placebo:Young-Drug   Old-Placebo:Old-Drug   Young-Placebo:Old-Drug   Old-Placebo:Young-Drug So I tried:               pH = glht(ANC, linfct = mcp(treatment*age="Tukey"))     Error: unexpected '=' in "ph = glht(ANC_nback, linfct = mcp(treat*age="      And:               pH = TukeyHSD(ANC)     Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep3'     In addition: Warning message:     In replications(paste("~", xx), data = mf) : non-factors ignored: baseline      Does anyone know how to resolve this? Many thanks!


title: Why is it necessary to use ML estimation instead of REML to compare multilevel linear models?

text: In Discovering Statistics Using SPSS 4e, Andy Field writes on p835 that: > SPSS gives you the choice of two methods for estimating the parameters in > the analysis: maximum likelihood (ML), which we have encountered before, and > restricted maximum likelihood (REML). The conventional wisdom seems to be > that ML produces more accurate estimates of fixed regression parameters, > whereas REML produces more accurate estimates of random variances (Twisk, > 2006). As such, the choice of estimation procedure depends on whether your > hypotheses are focused on the fixed regression parameters or on estimating > variances of the random effects. However, in many situations the choice of > ML or REML will make only a small difference to parameter estimates. Also, > if you want to compare models you must use ML. I'm most interested in the last sentence. Why is it that to compare models you must use ML and not REML?


title: curse of dimensionality & nonparametric techniques

text: I have seen it many times in a number of articles that nonparametric techniques are subject to the curse of dimensionality, which may lead to the failure of these methods. Why does this happen? Could give me an example? Thanks!


title: Statistics on mathematical "data"?

text: Lets say there is a theoretical relationship you want to prove over all values of some variable. For example, $F(x)$ increases with $x$. You are unable to come up with a general theoretical proof. However, you can calculate $F(x)$ for any specific value of $x$, and the expected relationship holds. Presumably, you can then run some sort of regression on the "data" provided by a whole bunch of the values you specifically calculated. (using $x$ as the IV and $F(x)$ as the DV). Then, you can show that an increasing relationship is supported by the "data" and try to come up with a best fit, etc. This is sort of like "comparative statics" only taking it one step further -- trying to prove the direction (even shape) of a relationship, rather than just the ordinal relationship of two values. Is there a name for this type of process -- running statistics on mathematical solutions, vs. measured data? Is this a convincing way to argue a theoretical point? Isasmuch as running an experiment and collecting data values and running statistics on that is? The final goal is not a theoretical contribution on $F$'s movement with $x$, but rather a model that provides a useful takeaway (e.g. since $F(x)$ increases with $x$, it is advisable to try to minimize $x$).


title: What is ridge regression?

text: I just need a simple explanation of what exactly ridge regression is so I can have a decent intuitive understanding of it. I understand it's about applying some sort of penalty to the regression coefficients... but beyond that I'm a little confused about how it is different from other kinds of regression which implement penalties. In what case should you use ridge regression as opposed to some other kind of regression?


title: Why does the event $\{X \leqslant t\}$ is equal to the sample space when $t \geqslant c$ and $X=c$?

text: This is what is written in my notes. > Constants are RVs: for $X=c=const$, one has $\\{X \leqslant t\\} = > \begin{cases} \emptyset, t<c \\\ \Omega, t \geqslant c\end{cases}$ Why does the event $\\{X \leqslant t\\}$ is equal to the sample space when $t \geqslant c$ and $X=c$?


title: Why are the Borel subsets on $\mathbb R$ a $\sigma$-algebra?

text: I am a newbie in mathematical statistics and haven't learnt any group theory before. My lecture notes are too brief. How can the Borel subsets on $\mathbb R$ satisfy A.2 and A.3 of a $\sigma$-algebra ? A.2 is closure under complement. A.3 is under countable unions.


title: Confused about comparison of Pearson coefficients

text: I am slightly confused on a paper I am reading. They have a table of calculated Pearson coefficients as shown here: $$ \begin{array} . & I & II & III\\\ A + B & 0.8 & 0.65 & 0.90\\\ A + C & 0.73 & 0.46 & 0.83\\\ B + C & 0.61 & 0.76 & 0.68\\\ \end{array} $$ They claim that the weakest correlation recorded was between $B+C, \,\, I$ and $B+C, \,\, III$. $A,B,C$ are different imaging modalities. Are they stating that these are the two closest values in one comparison (i.e $B+C$)? I am confused because they are now correlating correlations! I thought the the weakest correlations would be $B+C, \,\, I$ and $A+C, \,\, II$, as they have the smallest values, but perhaps I did not understand their statement. The paper can be found here: http://online.medphys.org/resource/1/mphya6/v40/i4/p041707_s1?bypassSSO=1 Page 4, section III.B


title: Naive Bayes Python implementation differences

text: Am currently using Naive Bayes for a multi labelled text document classification problem. But I would like to know the differences (advantages and disadvantages) of using SkiLearn Naive Bayes or NLTK? SkiLearn seems to have flexibility in parameter setting more than NLTK, but I'm think my thoughts on this being the only difference is quite naive. Also what cases would you select one over the other? Any help appreciated.


title: Feasible Generalized Least Square in R

text: I am studying the factors influencing the annual salary for employees at a undisclosed bank. The regression model that I have decided to employ is as follows: \begin{equation} Y_{k}=\beta_{1}+\beta_{2}E_{k}+\beta_{3}D_{gk}+\beta_{4}D_{mk}+\beta_{5}D_{2k}+\beta_{6}D_{3k}+\varepsilon_{k} \end{equation} where $Y_{k}$ is the logarithm of annual salary, $E$ is the number of years of education, $D_{g}$ is a gender dummy, $D_{m}$ is a minority dummy, and where \begin{equation} D_{2}=\begin{cases} 1 &\text{Custodial job} \\\ 0 & \text{Otherwise} \end{cases} \end{equation}   and \begin{equation} D_{3}=\begin{cases} 1 &\text{Management job} \\\ 0 & \text{Otherwise} \end{cases} \end{equation} As you know, whenever one deals with GLS, $\Omega$ will almost surely be unknown and thus have to be estimated. In general there are $\frac{n(n+1)}{2}$ parameters to be estimated, which makes it pretty impossible to come up with a viable estimation out of $n$ observations. This is usually counteracted by imposing some structure on $\Omega$. In my case, I would like to make the assumption that the disturbance terms $\varepsilon_{k}$ in the above regression model have variance $\sigma_{i}^{2}$ for $i=1,2,3$, according to whether the $i$-th employee has a job in category 1,2, or 3 respectively. Now, we may introduce the transformations $\gamma_{1}=\log (\sigma_{1}^{2}),\gamma_{2}=\log(\sigma_{2}^{2}/\sigma_{1}^{2})$, and $\gamma_{3}=\log(\sigma_{3}^{2}/\sigma_{1}^{2})$ so as to enable us to formulate the following model for \begin{equation} \sigma_{k}^{2}= \exp \\{ \gamma_{1}+\gamma_{2}D_{2k}+\gamma_{3}D_{3k} \\} \end{equation} Since $\hat{\beta}_\rm{OLS}$ is a consistent estimate of $\beta$, even under the assumption of heteroscedasticity, we have that $\hat{\beta}_{\rm OLS} \xrightarrow[]{p}\beta$ as the number of observations increase. We may therefore argue that $e_{k}^{2} \approx \sigma_{k}^{2}$, and so we can regress upon information that we already possess. **Summary of procedure** **(1)** Calculate the OLS estimate. **(2)** Calculate the OLS residual $\textbf{e}=\textbf{Y}-\textbf{X}\hat{\beta}$ **(3)** Calculate the OLS estimate of $\gamma$ from $e_{k}^{2}=f_{\gamma}(Z_{k})+\overline{\varepsilon}_{k}$. **(4)** Calculate the FGLS estimate as the GLS estimate with $\hat{\Omega}=\Omega(\hat{\gamma})$ in place of $\Omega$. > What I would like to know is whether or not one can perform this estimation > using a known function in R, say `gls`? If the answer is yes, then how > exactly should I write to ensure that that my heteroscedasticity assumption > is taken into account? Thanks for taking the time! Have a great day!


title: Multiple Regression

text: I need help in this area for my thesis. I have 2 time lines in my research in which different variables are measured. The only variable that is the same for both timelines is on stress. It was suggested to me that I can use multiple regression on a variable from time 1 and another variable in time 2 (as IVs) and testing it with stress just using time 2 (DV). Is that valid? Do I have to ensure any proceedings before I conduct the test? Or can I head on to check on assumptions of the data and I can run the test already? Do advice as Time 1 response is different from Time 2 response with missing data from drop outs. Thank you in advance!


title: How to calculate quartiles with only standard deviation and mean assuming normal distribution?

text: I am trying to understand how, using only the standard deviation and mean, you are able to determine the first and third quartiles on a normal distribution. I get that the area under the curve equals one; I understand that Q1 accounts for 25% of the area, and Q3 75%; but I do not understand how to leverage this information to calculate (with my TI-89) where Q1 and Q3 are... I was able to reasonably estimate these values by looking at a z-score table for areas representing 25% and 75%, and I got z-scores of ±0.68. However, I need a greater degree of accuracy than manually estimating from a table, but I do not know how to go about this formulaically or through my TI-89.


title: 1-sided two sample t-test with covariate adjustment

text: I have a biological problem which involves using a 1-sided two-sample t-test. To simplify the problem, let’s assume an experiment comparing quantities in a control (C) and a treatment (T) group, each having 3 replicates, called C1, C2, C3, T1, T2, T3. Here, **C1 and T1 are samples from the same subject, C2 and T2 are samples from the same subject, and so on.** Now I obtain a number for each of the following combinations. "within" = within group (i.e. both are C or both are T). Also, the number after C or T indicate whether they're from the same subject (sub) or not.               C1_C2 (within: different sub)     C1_C3 (within: different sub)     C1_T1 (between: same sub)     C1_T2 (between: different sub)     C1_T3 (between: different sub)     C2_C3 (within: different sub)     C2_T1 (between: different sub)     C2_T2 (between: same sub)     C2_T3 (between: different sub)     … (a total of 15 pairwise quantities)      We are essentially interested in comparing the quantities "within" group and the quantities "between" group, but **hoping to incorporating the "subject" effect somehow in the t-test**. Is there a way to distinguish “different sub” and “same sub” in the test? Thank you so much for your suggestions!!


title: zero inflated binomial data

text: I am working with presence/absence data that contains lots of zeros. What is the best method to model this. The only suggestions I can find for zero inflated data refer to count data not binomial/proportion data.


title: Can the normalizing constant of a PDF be infinite?

text: Given a probability density function (pdf), $$P(x) = \frac{f(x)}{\Sigma_x f(x)}.$$ Is this still a valid pdf if the normalizing constant in the denominator, $\Sigma_x f(x) = \infty$?


title: Effective validity of AUROC as performance measure: what about "very high" AUROC values?

text: The Area Under ROC curve (AUROC) is a quantity used to quantify performance of classifiers. I am currently interested in the most basic drawbacks of using AUROC as unique performance measure. I got a negative feedback with AUROC in some recent analysis as I met a classification problem performed through a logistic regression in presence of complete separation; I presume that in this case one simply obtains a perfect AUROC=1. Is this correct? In general, extremely high values of AUROC (I mean values $\geq 0.998$ ) would make me wonder whether there is something "wrong" or "artifical" with the model and/or data, as first feeling. It would be nice if you could confirm this feeling of mine and support it with examples. At this stage, the opposite thesis is also interesting to me.


title: Statistically comparing two different cache hit algorithms

text: I am comparing two different cache hit algorithms using four different training and testing sets. The training sets are related: 1 is a subset of 2 is a subset of 3 is a subset of 4. The same goes for the test sets. The result is two different matrices with cache hit ratios between 0 and 1, where higher values are better: Algorithm A                                 test set 1:  test set 2:  test set 3:  test set 4:     training set 1:   0..1         0..1         0..1         0..1     training set 2:   0..1         0..1         0..1         0..1     training set 3:   0..1         0..1         0..1         0..1     training set 4:   0..1         0..1         0..1         0..1      Algorithm B:                                 test set 1:  test set 2:  test set 3:  test set 4:     training set 1:   0..1         0..1         0..1         0..1     training set 2:   0..1         0..1         0..1         0..1     training set 3:   0..1         0..1         0..1         0..1     training set 4:   0..1         0..1         0..1         0..1      I now want to determine the best algorithm. How can this be done? I have limited understanding of statistical concepts. I don't know the distribution of the data. I thought about using a dependent sample t-test for each pair of cells in A and B, so that B is seen as the "treatment" and we then decide if the "treatment" was benefitial or not. H0: There is no significant difference between the algorithms   H1: There is a significant difference between the algorithms Is this method even remotely valid? Any help much appreciated!


title: Finding the mean difference of two data sets in R, within a confidence interval

text: I'm analysing benchmark data for some of my R functions, and I'm trying to compute a statistic describing how many times faster my function - group_one - is than a control function - group_two - as a multiplier (group_one is $\alpha$ times faster than group_two). The data are assumed to be normally distributed, of equal length. The variance is assumed to be different between groups.               test_data <- list(         group_one = c(9387298, 9447571, 10343495, 10232936, 10296492, 10273234, 16732716, 10527317, 10675171, 9434651)         group_two = c(23721104, 39465146, 31087543, 39851996, 38027882, 44277978, 32448261, 39950612, 47677503, 38887976)          # how many times smaller is the mean of the first group than the second?     multiplier <- mean(test_data$group_one) / mean(test_data$group_two)      The variance of the the benchmark data is typically about 15% of the mean, so I'm trying to find a 95% confidence interval in which the multiplier lies (something of the form "group_one is 10 - 12 times faster than group_two") Does anyone know an R function (preferably in the base library) that could be used to compute this interval? @COOLSerdash: I tried t.test, but it doesn't directly give an confidence interval in which the population mean lies, and I don't know how to infer that interval based on the t test. @Peter Flom: I added more realistic data above; the initial "data" were just random numbers. I can see that the ratio of the two distributions might be too sensitive to usefully gauge how much faster group_one is compared to group_two; if an alternative statistic would be more useful I'm open to suggestions.


title: Linear SVM C optimisation data, how to partition into train, model construction, test

text: My aim is to find the best C for a linear SVM classification using libsvm, I have 120 instances in total and 2 classes which I want to classify. I have question regarding partitioning a the dataset into these parts: 1) training for optimisation during parameter C search 2) training for model construction 3) testing (representing future incoming data) Using cross-validation and grid-search, I'd like to find the optimal C for my dataset. There seem to be different suggestions on how to partition data into these three parts I could first divide my data into training (50%) and testing (50%). From the training, I use one half (25% of total data) for cross-validation for optimisation of C. Once I found the best C, I use THE OTHER HALF of the training data (25% of total data) , to construct the model using the best C. I then validate the model against the 50% testing data, the resulting accuracy is the final accuracy. This way, neither the data that went in for the model- construction nor the testing data have ever been used during optimisation. Use 50% for training and 50% for testing. The entire 50% for training is used for cross-validation for optimisation of the parameter C. Once I found the best C, I train on the full 50% of the training data to obtain the model, and validate it on the test data. I think the advantage of 1. is that neither model construction and testing data were in the optimisation process, hence there is less of an overfitting problem. I read differing accounts online though, so any input would be greatly appreciated. Does the method depend on how much data there is available? In my case, I don't have that many instances, so I am not sure if method 1. is recommended. Reference to other relevant sources/textbooks would also be greatly appreciated.


title: Correcting for multiple comparisons on planned contrasts with an ANCOVA

text: I have done an experiment with 7 treatment conditions designed to reduce alcohol intake, and a control condition. The dependent variable is alcohol intake after treatment, and I have run an ANCOVA with baseline alcohol intake as a covariate. I found no effect from the omnibus ANCOVA F-test, but wanted to investigate specific comparisons between the treatment conditions (which was the main aim of the research). I therefore used planned contrasts to look at two specific questions:   * The average of all treatment conditions vs. control condition = no effect.   * Treatment 1 vs. all the other treatments individually. e.g. treatment 1 vs. treatment 2, treatment 1 vs treatment 3, treatment 1 vs treatment 4 etc. I have found a significant difference between treatment 1 and treatment 6 (which I had hypothesised - but not based on _strong_ evidence, more of a hunch) However, I have not applied any correction for an inflated type 1 error rate for doing these multiple comparisons. Everything I've read on this topic so far has been pretty vague as to how to use corrections, e.g. "oh its ok as long as you're only doing 'a few' planned contrasts or if you have an a priori reason for it" - but thats not very helpful, how many is 'a few' ? My question is: Do you think I need to apply an adjustment for multiple comparisons (I have done 7 in total) and if so, which adjustment? I've read that you can only use certain adjustments when using an ANCOVA (as opposed to an ANOVA). Dunnett's test looks like it could be right, but I don't know if it is possible to use Dunnett's test with an ANCOVA? Any help (or even pointing me in the direction of some more useful reading) would be much appreciated.


title: Variance of powers of a random variable

text: Is it possible to derive a formula for variance of powers of a random variable in terms of expected value and variance of X? $$\operatorname{var}(X^n)= \,?$$ and $$E(X^n)=\,?$$


title: How can I create binary variables from Likert scale for binary logistic regression model?

text: I used 5-point Likert-type scales as a dependent variables (1 strongly agree, 2 agree, 3 neutral, 4 disagree, 5 strongly disagree), but, as you are aware, using binary logistic regression requires me to create a binary variable from these 5-level responses. So, strongly agree and agree would become category 1, and disagree and strongly disagree would become category 2, while we exclude the neutral response option (although this significantly cuts down my sample size). One more issue in this regard, my questionnaire has 12 variables and each variable is measured by a group of items. For example, I have education (measured by 7 items), Gender inequality (measured by 6 items), trust (measured by 6 items), career development (measured by 7 items) and so on. Now, and in order to create the binary category do I need to take the high and the low score of each variables and find the mean and standard deviation, and exclude the cases in the middle? but this again will significantly cuts down my sample size? Finally is it better to use multinomial logistic regression instead of binary logistic regression when dependent variables are measured by 5-point Likert scale, which means there are 5 groups with each variable?


title: Model building based on preliminary data analysis

text: I have a extremely simple data set which consists of one independent variable (five sampling sites as factors) and the outcome variable which is disease status (positive/negative) in crabs from each site. These sites are located at varying distances from a sewage outfall point. My hypotheses are that 1) prevelence is different among sites, and 2) disease prevalence varies as a function of distance from an identified point-source pollution area. Hypothesis 1 was tested using a chi-square test for trend. Regarding hypothesis 2, I built a logistic regression model. Instead of using as a reference group the site located farthest away from the pollution source (least exposed crabs), I resorted to a preliminary analysis and used the site with the lowest prevalence as the reference group. Having said that, I did not choose my baseline category (against all other are compared) _a priori_ , but only after knowing that crabs from that site are the least affected ones. It made biological sense to me to build the model that way; however, I am not sure that I am not violating any assumption here by not including as reference category the one which I thought would have the lowest disease prevalence in the first place. Am I making a big mistake here? Thank you in advance.


title: What is the name of the theory (if there is one) that states lottery winners are more likely to tell others about their lottery entry?

text: I've previously read about a theory that I remember (correctly or not) being called the "Winning Lottery Theory" which is essentially the following: > An individual hears about disproportionately more occurrences of other > individuals entering a lottery and winning, than they do of others entering > a lottery and not winning. This gives the individual a distorted perspective > of the odds of actually winning. The same principle applies to people's testaments of products or services that work for them as opposed to not working for them. The opposite effect also applies, if a detrimental effect occurs. It seems particularly applicable to subjects likely vaccination (eg causing autism), and pseudoscience (eg my back stopped hurting after my homeopathy drops) etc. **Anyway, the question is this: is there a more standard name for this effect?** Googling "Winning lottery theory" came up with very little. The closest I could find was this article mentioning a "Lottery-Winner Theory", which returns a few more Google results but still not many.


title: Getting the MVUE of $\sigma^4$ for normal distribution

text: If x$_{1}$, x$_{2}$,...,x$_{n}$ ~ N($\theta$, $\sigma$$^{2}$), based on the fact that Y$_{1}$ = $\sum$(x$_{i}$) and Y$_{2}$ =$\sum$(x$_{i}$)$^{2}$ are complete and sufficient statistics, it can be shown that (not sure how) MVUE of $\sigma$$^{2}$ = $\sum$(x$_{i}$- $\bar x$ )$^{2}$ / (n-1) If I need to find the MVUE of $\sigma$$^{4}$, can I say it's simply [$\sum$(x$_{i}$- $\bar x$ )$^{2}$ / (n-1)]$^{2}$ ? I think the answer is no because E[A] $.$ E[A] does not neccessarily equal E [A$^{2}$]. Can someone help me confirm please? Thanks


title: What is the meaning of the term "enrichment" when performing cross-validation?

text: Trying to understand a discussion of a 5-fold cross-validation process to validate a predictive model and its results, there is a particular phrase which has me stumped, i.e.: > The predictions of the algorithm are significantly enriched with the left- > out test data. **What I think it means**. Trying to understand this, the way which seems more natural for me to interpret the statement is: "When we add the left-out test data, the model is better," which seems pretty obvious. To me, this would imply that the model needs that data - and _failed validation_. **What it seems to me it should mean**. However, in the context of validation, I would expect something like: "There is a high degree of similarity between what our model predicted without the left-out data and said left-out data". To me, this would imply that the model _is good at prediction_. I believe I am not understanding the concept of enrichment. After looking it up, the meaning is, in broad terms, "improvement". If this is so, and therefore my first interpretation is correct, why is this important for validation? What does "enrichment" mean in this context? Thanks!


title: Why do I get clipping when using tanh activation function?

text: I implemented a rather simple MLP NN as a part of my project. I'm just testing it on a sinusoid right now and I expect this network's output to follow the sinusoid without issues. Some network details: 32 input nodes 10 hidden nodes I use tanh activation functions in both layers I trained it using both BP and Prediction Error Algorithm and I see the same issue - I get a clipping of the tip of the function. Here is an image to illustrate: ![Clipped image](http://i.stack.imgur.com/aHy7O.png) If I make the target sinusoid with an amplitude of 0.9, it works fine: ![Good prediction image](http://i.stack.imgur.com/xz7sS.png) I have tried to crank up the learning rate up to 0.5, but it doesn't help (it just looks noisier). I think it's something I'm missing in my derivatives perhaps. Here are my update equations (please ignore absence of transposes): $$ y = tanh(H*tanh(Wx)); $$ $$ e = y - y_d $$ $$ dH = \mu_1 * e * (1-tanh(H*tanh(Wx)^2)) * tanh(Wx) $$ $$ dW = \mu_2 * x * e * H * (1-tanh(Wx)^2) * (1-tanh(H*tanh(Wx)^2)) $$ $$ H = H - dH $$ $$ W = W - dW $$ I know it does the clip because it doesn't reach a value of tanh that would output 1, (because tanh is nonlinear). Maybe I should squish my activation function to make it narrower? That's it... Any thoughts would be appreciated.


title: Using Weighted Least Squares with Robust Standard Errors

text: Weighted least squares (WLS) and robust standard errors are sometimes presented as alternative approaches for obtaining reliable standard errors of estimates of regression coefficients in the presence of heteroscedasticity. However, I notice that my software (gretl) offers robust standard errors as an option when using WLS. A situation in which it seems this might be useful is where, in a regression of Y on X, there is a clear reason for heteroscedasticity, for example a scale effect such that larger values of Y are expected to be associated with larger variances. One might then use WLS, giving a higher weighting to observations with smaller Y (or, perhaps better, to observations with smaller E[Y | X], as inferred from an initial OLS regression). However, it might be found that the WLS residuals suggested some remaining heteroscedasticity that the weighting had not eliminated. This would suggest that the standard errors estimated by WLS might not be entirely reliable, and to address this one might opt for robust standard errors (rather than attempting to do so via some more complex weighting pattern). Question: Assuming the number of observations is reasonably large (over 100, say), are there any pitfalls in using WLS with robust standard errors when estimating standard errors of regression coefficients?


title: Urn with non-uniform probabilities

text: An urn contains N-1 red and 1 green ball. Each ball has an associated weight. If each ball is drawn (without replacement) with a probability proportional to how much its weight contributes to the urn, what is the expected number of attempts required to get the green ball? Example: 2 red balls of weight 0.3 and 0.4, and green ball of weight 0.3. In first attempt, `Pr(red1)=0.3`, `Pr(red2)=0.4`, `Pr(green)=0.3`. Say, the red ball with weight 0.3 is chosen in the first attempt.   Then, in the next attempt, `Pr(red2)=0.4/(0.4+0.3)` and `Pr(green)=0.3/(0.4+0.3)`. It becomes relatively difficult to keep track of the probabilities if there are more red balls.


title: Comparison of entropy and distribution of bytes in compressed/encrypted data

text: I have some question which occupies myself for a while. The entropy test is often used to identify encrypted data. The entropy reaches its maximum when the bytes of the analyzed data are distributed uniformely. The entropy test identifies encrypted data, because this data has a uniform distribution - like compressed data, which is classified as encrypted when using the entropy test. Example: The entropy of some JPG file is 7,9961532 Bits/Byte, the entropy of some TrueCrypt-container is 7,9998857. This means with the entropy test I cannot detect a difference between encrypted and compressed data. BUT: as u may see on the first picture, obviously the bytes of the JPG-file are not distributed uniformely (at least not as uniform as the bytes from the truecrypt-container). Another test can be the frequency analysis. The distribution of each byte is measured and e.g. a chi-square test is performed to compare the distribution with a hypothetic distribution. as a result, I get a p-value. when i perform this test on JPG and TrueCrypt-data, the result is different. The p-Value of the JPG file is 0, which means that the distribution from a statistical view is not uniform. The p-Value of the TrueCrypt-file is 0,95, which means that the distribution is almost perfectly uniform. My question now: Can somebody tell me why the entropy test produces false positives like this? Is it the scale of the unit, in which the information content is expressed (bits per byte)? Is e.g. the p-value a much better "unit", because of a finer scale? Thank you guys very much for any answer/ideas! JPG-Image ![enter image description here](http://i.stack.imgur.com/4LE3G.png) TrueCrypt-Container ![enter image description here](http://i.stack.imgur.com/8gxfF.png)


title: ANOVA claims that model is singular

text: I have a problem with a 3 way ANOVA with repeated measures. The design is the following: "localization" is the dependent variable, than I have the following within subjects factors "material" (4 levels), "diffusion technique" (6 levels) and "typology" (2 levels). This is the formula I used:               aov_Material = aov(abs(Localization) ~ Diffusion_Type*Material*Typology +                        Error(Subject/(Diffusion_Type*Material*Typology)), data=scrd)     summary(aov_Material)      And this is the output I get as result:               Warning message:     In aov(abs(Localization) ~ Diffusion_Type * Material * Typology +  :       Error() model is singular     > summary(aov_Material)          Error: Subject               Df Sum Sq Mean Sq F value Pr(>F)     Residuals 11   2117   192.5                         Error: Subject:Diffusion_Type                    Df Sum Sq Mean Sq F value Pr(>F)     Diffusion_Type  5   20.2   4.032   0.423   0.83     Residuals      55  523.7   9.522                         Error: Subject:Material               Df Sum Sq Mean Sq F value   Pr(>F)         Material   3   1967   655.6    21.4 7.03e-08 ***     Residuals 33   1011    30.6                          ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Error: Subject:Diffusion_Type:Material                              Df Sum Sq Mean Sq F value Pr(>F)     Diffusion_Type:Material  15  170.1  11.338   1.308  0.202     Residuals               165 1430.3   8.668                         Error: Within                Df Sum Sq Mean Sq F value Pr(>F)     Residuals 288   2888   10.03         My questions are:   1. Why there is the warning (and, BTW, what does it means that the model is singular)?   2. Why there is not the effect for typology listed in the results? How can I understand if there is a difference between the two levels of typology?


title: Interpreting p-value significance from Pearson correlation

text: I used matlab corr() function to identify correlation of 236 samples. Pearson correlation is selected, and the output return r and p-value. Two sets of samples returned different r & p-value. May I know how to interpret the significance of correlation with the below results? (a) The data has strong negative correlation, and it's significant as p-value is a lot lesser than 0.05 ( p << 0.05 )               r = -0.9383        p = 6.7415e-110      (b) the data has weak positive correlation, and it's insignificant as p-value > 0.05.               r = 0.06800       p = 0.2981        Am I right?


title: Question on leave one out and stratified 10-fold cross validation

text: I am confused with the answers to the questions below Assume that we have a dataset D with 100 examples, 50 of which belong to the class ’good’ and 50 belong to the class ‘poor’. Assume further that we have a very naïve learning algorithm (L) that produces a model that simply predicts the majority class of the training examples (if there is no majority, it always predicts ‘good’) 1) Performing stratified 10-fold cross validation on D will give the same result as when using half the examples of each class for training and the other halves for testing. Answer> True 2) Leave-one-out cross validation will give a better accuracy estimate for how L will perform on examples that are independently sampled from the same distribution as D than when using half the examples for training and half the examples for testing.Answer> False. May I know why the first answer i true and the second is false?


title: Diagnosing why MICE is crashing R when attempting to impute multilevel data

text: I have never had problems with R crashing before. I am using the `mice` package (mice 2.13) to perform multiple imputations. The code works fine on some subsets of the data, but when I run it on other subsets, R crashes (not immediately - after some time). From the output in R just before it crashes, I believe it is using the `2l.pan` method of imputation (from the `pan` package) I have run update.packages() already. How can I diagnose this problem ?               Problem signature:       Problem Event Name:   APPCRASH       Application Name: Rgui.exe       Application Version:  2.151.59607.0       Application Timestamp:    4fe47a63       Fault Module Name:    R.dll       Fault Module Version: 2.151.59607.0       Fault Module Timestamp:   4fe47a4e       Exception Code:   c0000005       Exception Offset: 0000000000032ec8       OS Version:   6.1.7601.2.1.0.256.4       Locale ID:    2057       Additional Information 1: 7782       Additional Information 2: 77823beb5887f451c3dd7ae4fe931995       Additional Information 3: 4491       Additional Information 4: 4491b41bf90894717964f5eef2cccd84      **Update** I have managed to create a reproducible example, with data:               require(foreign)     require(mice)     require(pan)          dt.fail <- read.csv("http://goo.gl/pg8um")     dt.fail$X <- NULL          dt.fail$out <- as.factor(dt.fail$out )     dt.fail$grp<- as.factor(dt.fail$grp)     dt.fail$v1<- as.factor(dt.fail$v1)     dt.fail$v2<- as.factor(dt.fail$v2)     dt.fail$v3 <- as.factor(dt.fail$v3)     dt.fail$v7<- as.factor(dt.fail$v7)     dt.fail$v8 <- as.factor(dt.fail$v8)     dt.fail$v9 <- as.factor(dt.fail$v9)     dt.fail$v11 <- as.factor(dt.fail$v11)     dt.fail$v12 <- as.factor(dt.fail$v12)               PredMatrix <- quickpred(dt.fail)     PredMatrix['CTP',] <- c(1,-2,0,0,0,0,0,0,0,0,1,0,1,1,0,2)                    impute = mice(     data=dt.fail,          m = 1,          maxit = 1,         imputationMethod = c(         "logreg",   # out         "",      # grp   ----> cluster grouping factor         "pmm",  # v1         "polyreg",  # v2         "logreg",   # v3         "pmm",  # v4         "logreg",   # v5         "logreg",   # v6         "polyreg",  # v7 ----> auxilliary         "polyreg",  # v8 ----> auxilliary         "polyreg",  # v9 ----> auxilliary         "polyreg",  # v10 ----> auxilliary         "",     # v11 ----> complete         "",     # v12 ----> complete         "2l.pan",   # CTP ----> multilevel imputation         ""),        # const ----> needed for multilevel impuitation          predictorMatrix = PredMatrix, seed = 101     )      And for completeness, here is the predictor matrix I was using:                   .     out grp v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 CTP const     out     0   0  0  0  0  0  1  1  0  0  0   0   0   0   0     0     grp     0   0  0  0  0  0  0  0  0  0  0   0   0   1   0     0     v1      0   0  0  0  0  0  0  0  0  0  0   0   1   1   0     0     v2      0   0  0  0  0  1  1  1  0  1  0   0   1   1   1     0     v3      0   0  0  0  0  1  1  1  0  1  1   0   1   1   1     0     v4      0   0  0  1  1  0  1  1  0  1  1   0   1   1   1     0     v5      1   1  0  0  0  0  0  1  0  1  0   0   1   0   0     0     v6      1   1  0  1  0  1  1  0  0  1  0   0   1   0   0     0     v7      0   0  0  0  0  0  1  1  0  1  0   0   0   1   0     0     v8      0   0  0  0  0  0  1  1  0  0  0   0   1   1   0     0     v9      0   0  0  0  1  1  1  1  0  1  0   0   1   1   1     0     v10     0   0  0  0  0  0  1  1  0  1  0   0   1   1   0     0     v11     0   0  0  0  0  0  0  0  0  0  0   0   0   0   0     0     v12     0   0  0  0  0  0  0  0  0  0  0   0   0   0   0     0     CTP     1  -2  0  0  0  0  0  0  0  0  1   0   1   1   0     2     const   0   0  0  0  0  0  0  0  0  0  0   0   0   0   0     0


title: How many samples are needed to estimate a p-dimensional covariance matrix?

text: In general, how many points are needed to estimate a p-dimensional covariance matrix? Does it depend on how the data are spread out across the different dimensions? Does it depend on the true distribution of the data? Thank you!


title: Different Mean Square partitions in an unbalanced bifactorial ANOVA (with random factor) between R and Statistica

text: I am trying to extract variance components for selection and chance in a bifactorial design with **Generation as a fixed factor** and **Replicate as a random term,** for early fecundity. Since I am using the individual level (because I want to make confidence intervals for Chance), the sample size is not equal between Generations so it's an unbalanced design and it's is not a repeated measures ANOVA because my individuals are different. The problem is that when I compare the results between Statistica and R the Mean Squares(MS) for the Generation and Replicate Factors are different (but the MS for Gen*Rep and the residuals are not!) My data consist of 2 terms:   * Generation - 2 level fixed factor (6 and 11)   * Replicate - 3 level random factor (Ad1, Ad2, Ad3) The model I am using is               g <- Anova (aov(Fec~ Gen+Rep+ Gen*Rep), random=~1|Gen*Rep)      I am not sure if this is how to phrase that I want the terms with Replicate to be random (maybe that's one of my errors when comparing with the Statistica results?).   * Might the differences be due to the unbalanced design? Do R and Statistica have different ways of handling the different number of individuals in each replicate and generation?   * Why aren't all MS equal? Is there a way to make the results be the same?    * How can I indicate that a factor (and all its interactions) are random?    * How can I tell R that I want it to test the Generation factor against the Gen*Rep term or to put it in another way to make it test the significance of the Gen factor against the heterogeneity between my replicates?


title: Determine accuracy of model which estimates probability of event

text: I'm modelling an event with two outcomes, a and b. I have created a model which estimates the probability that either a or b will happen (i.e. the model will calculate that a will happen with 40% chance and b will happen with 60% chance). I have a large record of outcomes of trials with the estimates from the model. I would like to quantify how accurate the model is using this data - is this possible, and if so how?


title: t-test what needs to be normally distributed?

text: I'm aware that the t-test needs 'normally distributed data'. But take the variable `y`. When it is plotted without being split by `group`, it isn't normally distributed:               set.seed(1)     y <- c(rnorm(1000, 1), rnorm(1000, 5))     group <- c(rep("A", 1000), rep("B", 1000))     df <- data.frame(y=y, group=group)     library(ggplot2)     ggplot(df, aes(y)) + geom_histogram()      ![enter image description here](http://i.stack.imgur.com/pWJqp.jpg) But when `y` is split by `group`, it is normally distributed:               ggplot(df, aes(y)) + geom_histogram() + facet_grid(~group)      ![enter image description here](http://i.stack.imgur.com/sgWN8.jpg) Can anyone clarify if a variable only needs to be normally distributed after being split by group?


title: How can I transform time series data so I can use simpler techniques for fault prediction?

text: I know this is primarily a statistics site, so if I am off-topic, please redirect me. I have a system with pumps that sometimes break and need to be replaced. I would like to be able to predict the failures, and thereby give early warning to the people replacing the pumps. I have historical data for the pump process, such as flow, pressure, liquid height etc. I have only a small amount of experience in using machine learning techniques to classify data - basically I have followed and done the exercises of Andrew Ng's machine learning course on coursera, as well as Andrew Conway's Statistics One, - and I have never used machine learning to classify time series. I am thinking of ways I can transform the my problem so that I can use my existing knowledge on it. With my limited knowledge, I will not get a very optimal prediction, but I hope to learn from this, and for this problem, any small improvement in prediction is useful, versus just waiting for the faults to occur. My proposed approach is to turn the time series into a normal classification problem. The input would be a summary of a time series window, with mean value, standard deviation, max values etc. for each type of data in the window. For the output, I am not sure what would work best. One approach is that the output would be a binary classification of whether the pump failed within a certain time period from the end of the window or not. Another is that the output would be the time left before the pump fails, so not a classification, but a regression (in the machine learning sense) instead. Do you think this approach is likely to yield results? Is it a question of "depends on the domain and historical data". Are there better transforms (of both input and output) that I haven't considered, or is fault prediction based on time series data so different from more standard fault prediction, that my time would be better spent reading up on machine learning with time series?


title: Predicting variance of heteroscedastic data

text: I am trying to do a regression on heteroscedastic data where I'm trying to _predict_ the error variances as well as the mean values in terms of a linear model. Something like this: $$\begin{align}\\\ y\left(x,t\right) &= \bar{y}\left(x,t\right)+\xi\left(x,t\right),\\\ \xi\left(x,t\right) &\sim N\left(0,\sigma\left(x,t\right)\right),\\\ \bar{y}\left(x,t\right) &= y_{0}+ax+bt,\\\ \sigma\left(x,t\right) &= \sigma_{0}+cx+dt. \end{align} $$ In words, the data consists of repeated measurements of $y(x,t)$ at various values of $x$ and $t$. I assume these measurements consist of a "true" mean value $\bar{y}(x,t)$ which is a linear function of $x$ and $t$, with additive Gaussian noise $\xi(x,t)$ whose standard deviation (or variance, I haven't decided) also depends linearly on $x,t$. (I could allow more complicated dependencies on $x$ and $t$ – there isn't a strong theoretical motivation for a linear form – but I'd rather not overcomplicate things at this stage.) I know the search term here is "heteroscedasticity," but all I've been able to find so far are discussions of how to reduce/remove it to better predict $\bar{y}$, but nothing in terms of trying to _predict_ $\sigma$ in terms of the independent variables. I would like to estimate $y_0, a, b, \sigma_0, c$ and $d$ with confidence intervals (or Bayesian equivalents), and if there is an easy way to do it in SPSS so much the better! What should I do? Thanks.


title: Treatment mean comparisons in a mixed model with random and fixed effects

text: I have run a mixed model with one random effect, one fixed effect and a second order covariate effect. The random effect models the individual experimental units in my study. The LSMeans calculated for each experimental unit from the full model are quite different from the effects calculated running an individual model for each experimental unit by itself using only the fixed effect and second order covariate in the model. Is this to be expected? If so, why are the full model effects so different from the effects calculated with one experimental unit at a time?


title: Separating Base and Promotional Volume

text: I am working on a project where I have to separate base and promotional volume from the sales data. I have sales data for the last 4 years at week level. How can I separate base and promotional volume from only this much data? Is it possible to do by regression or any time series approach? I'd appreciate any feedback. Thanks.


title: How to calculate log-normal parameters using the mean and std of the given distribution

text: I have the mean (u) and the standard deviation (sd) of a continuous distribution (X). How do I solve for the mean (u_log) and standard deviation (sd_log) of the log of that continuous distribution (log(X))? I am looking to draw random numbers from a log-normal distribution. The problem is, most statistical packages have a log-normal-random-number- generator function that takes the mean of log(X) and sd of log(X) as inputs. I only have the mean of X and std of X. How to I solve for mean(log(X)) and sd(log(X)) given those inputs.


title: What do these residual patterns tell me?

text: I am doing a regression analysis with multiple variables and comparing it to a one-variable null hypothesis. The goal is to see which model provides a better explanation. The topic is information diffusion, hence the names DIFfusion, INDustry, etc. Also I must say that DIF is many times zero (has a long tail) and all variables go between 0 and 1. Anyways, the hypothesis can be translated unto: H_0: DIF ~ REL H_1: DIF ~ REL+COM*REL+IND*REL H_2: DIF ~ REL+COM*REL+IND*REL+SIZE1+SIZE2 REL multiplies the other two variables as theoretically they are related and together they 'should' give a better prediction of DIF. Now, the tricky things for me are two. First, (and most importantly), although R values are over 0.6 -good enough for this topic- when I see the residuals they tend to follow a pattern that is similar for all three hypotheses. I really don't know why or how to address it. ![Residuals](http://i.stack.imgur.com/XRh2q.png) ![QQplot](http://i.stack.imgur.com/JXeCI.png) ![histogram](http://i.stack.imgur.com/ZCvJe.png) Second, I have the intuition that the dependent variable (DIF) behaves like an iceberg and sea-level. Meaning that the lower the sea level (represented by SIZE1-2), I could see more of the shape of the iceberg, where the shape of the iceberg is given by the other variables. Have you encountered situations like this? How would you model/test it? Any advice? Using R by the way.


title: Guideline to select the hyperparameters in Deep Learning

text: I'm looking for a paper that could help in giving a guideline on how to choose the hyperparameters of a deep architecture, like stacked auto-encoders or deep believe networks. There are a lot of hyperparameters and I'm very confused on how to choose them. Also using cross-validation is not an option since training really takes a lot of time!


What effect does data averaging have on the variogram? To be specific, please see a simple example:               #Simulate a pure random walk data set, call this *PROCESS A*     n <- 1000  #number of data points     t <- 1:n #time     y <- cumsum(rnorm(n)) #data points          # Averaging at every 2 lags, call this *PROCESS B*     t2 <- apply(matrix(t, nrow=n/2, byrow=TRUE), 1, mean)     y2 <- apply(matrix(y, nrow=n/2, byrow=TRUE), 1, mean)          # Compute and plot the empirical variogram     require(geoR)     var <- as.geodata(data.frame(coords = t(rbind(t, rep(1, n))), data.col=y))     vario <- variog(var1)      var2 <- as.geodata(data.frame(coords = t(rbind(t2, rep(1, n/2))), data.col=y2))     vario2 <- variog(var2)      plot(vario$u, vario$v, type='b', pch=16, cex=.7)     points(vario2$u, vario2$v, col=2, type='b', pch=16, cex=.7)      Theoretically $\gamma(h) = \frac{1}{2} Var(y(t) - y(t+h))$ and for a linear variogram of a simple process without the nugget effect this is $\frac{1}{2}\sigma^2 |h|$ where $\sigma^2$ is the variance of the underlying process. Thus, we expect (1) the two variogram are parallel and (2) the proportion of the nugget effect of PROCESS A to PROCESS B is 4/3. Compare the plots and the Math: (1) is confirmed, but (2) is not. Why? Is this a simple example of the 'change of support' problem in geostatistics?? Thanking you in advance.


title: What statistical tests would be effective for this design?

text: Here the basics: Group 1 is the control group, they are given a specific website to use. Group 2 is the experimental group, they are given a similar website, with adjusted hypertext properties. After they complete the required actions, they are asked to complete a survey, (Likert Scale) which attempts to measure their feeling on the security and convenience of the website. This is the condensed version of the question. Just want to check. The IV: Hypertext The DV: Security and Convenience (Measured through Likert Scale) But which statistical tests in SPSS can I use to indicate a significant difference in the control and experimental group? And added advice would also be great. New to stats, so any help would be appreciated.


title: Data Mining study and prediction of a Dataframe in R

text: I'm new in the Data Mining World. I have a Dataset of 19 variables(some of them categorical). It is about execution time of different aplicattions. I have something like this) but with thousands of rows.                  app process subprocess ex_time   place start_date start_time session     1   A    A.a       A.1       2   Spain 01/02/2014   17:02:14    1234     2   A    A.b       A.2       2   Spain 01/02/2014   17:02:16    1234     3   A    A.c       A.2       3  France 01/02/2014   12:15:24    8374     4   B    B.d       B.4       4 England 02/02/2014   17:13:56   94374     5   B    B.e       B.5       5 Germany 16/02/2014   10:11:14    3873     6   C    C.f       C.7       6  Poland 25/01/2014   11:33:15   83739     7   C    C.g       C.9       0  France  9/04/2014   12:02:14    9437      Each app has different kind of processes but each subprocess can be in a different process(but must be into the same app) For example subprocess A.2 is inside the process A.b and A.c(and must be inside A). ex_time is the execution time. Session is a number for every session of the user( it means that in for example session 1234, the user executed the first 2 rows of the sample). How can I study the date base? How can I predict the execution time? What should I use to to that? Decision trees, clustering, (I've tried with linear regression using the subprocesses ordered by apps)?, Time series, Neuronal network, association rules? Thank you!


title: Intraclass Correlation Coefficient Confusion

text: Many times in statistics we are using a statistic in different ways. I think Intraclass correlation coefficient is one of them... Intraclass Correlation Coefficient can be a A) measure of inter-rater reliability, but it can also be used as B) measure of dependency between scores (preliminary test of a multilevel analysis). So the question(s) is/are: Do we use the same method to measure it?? For instance, in GUI programs such as SPSS or JMP do we use the same method to measure it??? Can you please refer to a url or a paper that shows how to use it as (B)? What is the range (general rule of thumb) for ICC for B??? Are there any references for that??? Thank you very much in advance for your contribution. I am looking forward to your reply.


title: Sensors categorization

text: I have several groups of sensor, each group send the data with different probability. I want to classify these groups into 3 categories ( first category has no data to send , the second category has one sensor to send and the third category has more than one sensor sending data ) how can I write the probability distribution function for these categories


title: How find a distribution for a set of data?

text: Hi I have a set of data(message with different size) in excel like (over 5000 entries):               message Name               MessageSize     message1                   0.5M     message2                   10.2M     message3                   2.1M     ...     message n-1                40.52M     message n                  12.12M      I have successfully generated a frequency table and a bar chart with it. e.g. Different range of size have how many message occurred ?               0 - 5M 4000     5 - 10M 454     5 - 15M 234     ..     ...     45-50M  2      Below are some statistics for this set of data:               Mean            Median          Mode            Skewness    Kurtosis             2.085443743         0.884527206 9.762702942 5.067491934 37.95900794           Standard deviation    min             max     3.642138843       0.006651878    49.12574959      But when I tried to draw a normal distribution via excel NORMDIST() it's has turned very ugly because it's probably not a Normal distribution ! As mean, median & mod are totally different ! How could I find out what kind of distribution the data belongs to ? or model a statistic model for it ? # Thanks a lot. I think stats here is a great community, lots of useful answers in a short while just like stackoverflow :)


title: What is the difference between conditional and unconditional quantile regression?

text: The conditional quantile regression estimator by Koenker and Basset (1978) for the $\tau^{th}$ quantile is defined as $$\widehat{\beta}_{QR} = \min_{b} \sum^{n}_{i=1} \rho_\tau (y_i - X'_i b_\tau)$$ where $\rho_\tau = u_i\cdot (\tau - 1(u_i<0))$ is a re-weighting function (called "check"-function) of the residuals $u_i$. In a paper by Firpo et al. (2009) the authors state that conditional quantile regression does not give interesting effects. They say that conditional results cannot be generalized to the population (in OLS we can alwas got from conditional to unconditional via the law of iterated expectations but this is not available for quantiles). This is because the $\tau^{th}$ unconditional quantile $y_i$ might not be the same as the $\tau^{th}$ conditional quantile $y_i |X_i$. If I understand correctly, part of the problem is that which covariates are included in $X_i$ has an effect on the ranking variable $u_i$ because inclusion of covariates splits the error into observed and unobserved components. I just cannot quite understand why this causes problems. Here are my questions:   1) What makes conditional and unconditional quantile effects different from each other?   2) How can I interpret the coefficients of conditional quantile regressions?   3) Are conditional quantile regressions biased? References   Firpo, S. et al. (2009) "Unconditional Quantile Regressions", Econometrica, Vol. 77(3), pp. 953-973


title: Casting a multivariate linear model as a multiple regression

text: Is recasting a multivariate linear regression model as a multiple linear regression entirely equivalent? I'm not referring to simply running $t$ separate regressions. I have read this in a few places (Bayesian Data Analysis -- Gelman et al., and Multivariate Old School -- Marden) that a multivariate linear model can easily be _reparameterized_ as multiple regression. However, neither source elaborates on this at all. They essentially just mention it, then continue using the multivariate model. Mathematically, I'll write the multivariate version first, $$ \underset{n \times t}{\mathbf{Y}} = \underset{n \times k}{\mathbf{X}} \hspace{2mm}\underset{k \times t}{\mathbf{B}} + \underset{n \times t}{\mathbf{R}}, $$ where the bold variables are matrices with their sizes below them. As usual, $\mathbf{Y}$ is data, $\mathbf{X}$ is the design matrix, $\mathbf{R}$ are normally distributed residuals, and $\mathbf{B}$ is what we are interested in making inferences with. To reparameterize this as the familiar multiple linear regression, one simply rewrites the variables as: $$ \underset{nt \times 1}{\mathbf{y}} = \underset{nt \times nk}{\mathbf{D}} \hspace{2mm} \underset{nk \times 1}{\boldsymbol{\beta}} + \underset{nt \times 1}{\mathbf{r}}, $$ where the reparameterizations used are $\mathbf{y} = row(\mathbf{Y}) $, $\boldsymbol\beta = row(\mathbf{B})$, and $\mathbf{D} = \mathbf{X} \otimes \mathbf{I}_{n}$. $row()$ means that the rows of the matrix are arranged end to end into a long vector, and $\otimes$ is the kronecker, or outer, product. So, if this is so easy, why bother writing books on multivariate models, test statistics for them etc.? It is most effective to just transform the variables first and use common univariate techniques. I'm sure there is a good reason, I just am having a hard time thinking of one, at least in the case of a linear model. Are there situations with the multivariate linear model and normally distributed random errors where this reparameterization does not apply, or limits the possibilities of the analysis you can undertake? Sources I have seen this: Marden - Multivariate Statistics: Old School. See sections 5.3 - 5.5. The book is available free from: http://istics.net/stat/ Gelman et al. - Bayesian Data Analysis. I have the second edition, and in this version there is a small paragraph in Ch. 19 'Multivariate Regression Models' titled: "The equivalent univariate regression model" Basically, can you do everything with the equivalent linear univariate regression model that you could with the multivariate model? If so, why develop methods for multivariate linear models at all? What about with Bayesian approaches?


title: Two-way ANOVA on not completely normal residuals

text: I have 2 factors, one (Page type) of 2 levels and the other (Intensity) of 4 levels. I did a 2x4 Full Factorial Design with 10 replications, so I did 80 experiments. Then I wanted to do a two-way ANOVA on these data. Correct me if I'm wrong, but what I have to check is if the residuals are indipendent, homoschedastic and normally distributed. The indipendence and the homoschedasticity is verified by the "Residual by Predicted Plot" because there are no visibile trends or increasing spreads. The normality assumption is not verified (second plot), but it seems not a big deal because: > Furthermore, normality is the least important assumption of a linear model > (e.g., an ANOVA); the residuals may not need to be perfectly normal. > (http://stats.stackexchange.com/a/28635/36058) It seems like they are _almost normal_. So, are my thoughts correct? Can I correctly do a two-way ANOVA on these data? (Full ANOVA: http://i.imgur.com/7tR3Kv2.png) ![enter image description here](http://i.stack.imgur.com/B9lct.png)


title: Loss rate calculation in a Poisson process

text: I have a bit a naive question about Poisson process. Suppose I have a station that receives a number of packets $n(i)$ each time slot $i$. If the number of received messages is less than $M$ then the $n(i)$ packages will all be sent. If $n(i) > M$ then $n(i)-M$ packets will be lost and $M$ packets will be sent. Let assume packets arrive at a Poisson rate $\lambda$. What is the mean of loss rate (total lost packages / total number of received packets)? I formulated the problem in this way:   * Let $N(t)$ be the total number of packets at time $t$, $N(t) = \sum_{i=0}^{t} n(i)$   * The total number of lost packets is $e(t) = \sum_{i=0}^{t} (n(i)-M) I_{(n(i) >M)}$ The loss rate will then be $\eta(t) = \frac{e(t)}{N(t)}$. The mean loss rate will be: $$\lim_{t \to \infty} E[\eta(t)]= E\left[\frac{\sum_{i=0}^{t} (n(i)-M) I_{(n(i) >M)}}{\sum_{i=0}^{t} n(i)}\right]$$ But I don't know how to calculate it. Any leads would be appreciated.


title: Leave one out cross validation and set seed in R

text: I am having trouble to understand where and how many times to use set.seed(). I am performing a leave-one out crossvalidation on a set with 120 classified variables to predict another 300 of unknown class. Should I set the seed at the beginning of the code or at the beginning of each iteration of the leave- one-out cross validation in order to guarantee a run that can be repeated? Thanks for any help and explanations.


title: Finding correlation coefficient

text: if I have A and B with the following known variables: with $E[A]$, $E[B]$ , $\sigma_{A}$ , $\sigma_B$ and correlation coefficient: $\rho_{AB}$ (assign numbers if you like) Say: $C=0.6A+0.4B$ Then how can I find $\rho_{AC}$? I know how to find: $E[C]$ ,$Var[C]$, but not how to find $\rho_{AC}$... I hope someone can help me out..


title: Is there a general rule about max nr of variables to use in (generalized) linear model?

text: I am running a generalized linear model on a dataset with 19 individuals and have 4 variables of interest. There are furthermore a number of interactions that might be interesting to look at. I was wondering if there is a general rule of thumb (with reference please) about the nr of terms you can use in a model based on the sample size. Thanks.


title: How to analyze Max Diff via R or SPSS

text: I am interested in utilizing Maximum Difference Scaling (Max Diff) for a survey, but apart from the Sawtooth software, I was not able to find a good answer on how to best analyze Max Diff. What is the best way to represent the result in my dataset, and which steps are used to analyze the data? I would prefer to use R or SPSS for analysis.


title: Logistic Regression - Adding interactions makes Independent variable statistically insignificant

text: My name is Abhi & I am trying to better understand logistic regression by solving a few practice problem. I am using R and RStudio as the development environment **_Problem Statement_**   Given the age, sex and class(first,second,third) for each passenger can you predict if he survived or died when the titanic sank **_Attempt1_**   Simple Logistic Regression. Use age,sex and passenger class directly. The formula (in R) `Survived~Pclass+Sex+Age` This gives fairly decent results - accuracy of 79% and all of independent variables are statistically significant **_Attempt2_**   Add interactions between age,sex and passenger class. Accuracy has improved to 80% but age is no longer significant. Also none of the new terms(age-sex,age- class,sex-class,age-sex-class) are statistically significant. The formula (in R) `Survived~Pclass*Sex*Age`      Can some one explain why this is happening? I can accept that the new terms may not be significant but why is age no longer significant? Any help would be much appreciated


title: How to change column names in data frame in R?

text: names(mydat)[c(name)]<-c("newname")       From this, I know that the column/variable name "name" of the data frame mydat is replaced with "newname". My question is if, I want to do this by a loop so that I will have some thing like: newname1 newname2 newname3 newname4 and so on, how do I do it? This is what did and it did not work:               for(i in 1:4){     names(mydat)[c(name)]<-c("newname"i)     }      Is there a way to code this? many thanks to all who could be of help. Owusu Isaac


title: Two forms of test statistics in likelihood ratio tests

text: From Bickel and Doksum's Mathematical Statistics ![enter image description here](http://i.stack.imgur.com/pdDbe.png) ![enter image description here](http://i.stack.imgur.com/wCZNu.png) I was wondering what "$\Theta_0$ is of smaller dimension than $\Theta = \Theta_0 \cup \Theta_1$" means? Note: In the original text, "$\theta_0$" is a typo and should be $\Theta_0$. Why when "$p(x, \theta)$ is continuous in $\theta$ and $\Theta_0$ is of smaller dimension that $\Theta = \Theta_0 \cup \Theta_1$", then $L(x)= \lambda(x)$? Thanks and regards!