title: Can i do a paired t-test on a pre- and post questionnaire but without pairing the same student for the pre and post?

text: I really suck at this so hope you could help me ve a group of people attending a training course, in order to measure whether the training course is effective in improving attendees' knowledge. We did the following:               attendees will complete a pre-test (TRUE or FALSE only) before the training     attendees attending the training     attendees will complete a post-test (TRUE or FALSE only) after training      Both pre and post tests are having exactly the same questions, However the test is NOT tagged with barcode so i have NO way to know if a certain pre and post questionnaire belong to the same individual Each paper has been scored out of 20 ,, and the scores i have entered in SPSS for every student for both the pre and post questionnaire So for this can i just use the paired t test and will the interpretation still be valid or is there an alternative i can use to find the significance of the intervention and whether it was effective or not? or should i do the t-test to the sum of the scores for pre and post? Thanks alot


When do we calculate the population parameter instead of sample statistics? If there is a this kind of case which statistical tool should we use population parameter or sample statistics? A car company wants to know the average age of cars of their brand that are still on the road. How would you define the appropriate population? Will the car company calculate a population parameter or sample statistic? Why?


title: Noise correlations depending on distance

text: If one has noise being generated over time $t\in[a,b]$, and the correlations of the value of the noise at times $t_0$ and $t_1$ turn out to be distributed according to a density function that depends only on the distance $|t_0-t_1|$ between them, then what is the distribution of the noise? How do the correlations determine the global distribution? Does it, for example, have to be equivalent to some version of the Ising model necessarily? Thanks a lot!


title: Lies, Damn Lies and Statistics

text: Is there something about statistics that lends itself to this sort of saying, or is it just that people will say anything to support their case, and this includes citing irrelevant or incomplete statistics?


title: Expected value of a random variable differing from arithmetic mean

text: I am a student who is taking a random processes class. I have seen that expected value of a discrete random variable is equal to the arithmetic mean of the distribution provided the values it takes. Is it true for all random variables irrespective of the distribution? Is there a case or example where expected value differs from the arithmetic mean? Secondly I think it applies only for discrete random variables. I think for continuous random variables, the pdf is zero at particular points. So in that case can I say that expected value is not equal to the mean of random variable?


title: Recreating traditional null hypothesis testing with Bayesian methods

text: I am trying to recreate (in R) a frequentist hypothesis testing in Bayesian from, by calculating Bayes factors of the null (H0) and alternative (H1) models. The model is simply a simple linear regression that tries to detect a trend in global temp. data from 1995 to 2009 (here). Therefore, H0 is no trend (i.e. slope = 0), or similary, the H0 model is a linear model with only the intercept. So I calculated the `lm()` of both models to arrive at negative log likelihood values that are significantly different. The p-value for the H1 lm() model is 0.0877. I also calculated this in a Bayesian way by using MCMCpack, and I get negative log likelihood values that are **super duper uber** different. Log likelihood values of 13.7 and 4.3 are about a 10000 fold difference in their likelihood ratios (where >100 is considered to be "decisive"). The means and sds of the estimates are very similar, so why am I getting such different likelihood values? (particularly for the Bayesian H0 model) I feel like there is a gap in my understanding on marginal likelihoods, but I can't pinpoint the problem. Thanks               library(MCMCpack)          ## data: http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3gl.txt          head(hadcru, 2)     ##  Year      1      2      3      4      5      6      7      8      9     10     ## 1 1850 -0.691 -0.357 -0.816 -0.586 -0.385 -0.311 -0.237 -0.340 -0.510 -0.504     ## 2 1851 -0.345 -0.394 -0.503 -0.480 -0.391 -0.264 -0.279 -0.175 -0.211 -0.123     ##       11     12    Avg     ## 1 -0.259 -0.318 -0.443     ## 2 -0.141 -0.151 -0.288          hadcru.lm <- lm(Avg ~ 1 + Year, data = subset(hadcru, (Year <= 2009 & Year >= 1995)))     hadcru.lm.zero <- lm(Avg ~ 1, data = subset(hadcru, (Year <= 2009 & Year >= 1995)))          hadcru.mcmc <- MCMCregress(Avg ~ 1 + Year, data = subset(hadcru, (Year <= 2009 & Year >= 1995)), thin = 100, mcmc = 100000, b0 = c(-20, 0), B0 = c(.00001, .00001), marginal = "Laplace")     hadcru.mcmc.zero <- MCMCregress(Avg ~ 1, data = subset(hadcru, (Year <= 2009 & Year >= 1995)), thin = 100, mcmc = 100000, b0 = c(0), B0 = c(.00001), marginal = "Laplace")          -logLik(hadcru.lm)     ## 'log Lik.' -14.55338 (df=3)     -logLik(hadcru.lm.zero)     ## 'log Lik.' -12.80723 (df=2)          attr(hadcru.mcmc, "logmarglike")     ##           [,1]     ## [1,] -13.65188     attr(hadcru.mcmc.zero, "logmarglike")     ##           [,1]     ## [1,] -4.310564      ![alt text](http://www.skepticalscience.com/images/HadCRUT_1995_2009.gif)


title: decomposition of independent random variables

text: If X and Y are independent random variables such that X = A + B and Y = C + D? Are the pairs (A, C), (B, C), (A, D), (B, D) also independent? By this I mean whether A and C are independent, B and C are independent and so on.


title: Including a nested factor as random effect in a GLMM

text: Good afternoon, I'd like to ask for advice on including a nested factor as a random effect in a GLMM. I've read other threads in this forum, but still am not able to answer my question. Any help is greatly appreciated! I tested children on 12 problem solving/tool use tasks which can be divided into 2 groups according to the frequency of the corresponding tool use behavior observed in wild great apes (low, high). Each child got a randomized set of 4 games, 2 low and 2 high frequency tasks. These are the variables: Dependent variable: Success (yes, no) Fixed effects: age (in months) sex (male, female) frequency (low, high) Random effects: subject tool use skill (rated by nursery staff, 1=below average, 2=average, 3=above average) task position (1-4) task (1-12) My main research question is whether the frequency variable is able to predict success. Task is nested within frequency, since task 1-6 = low and 7-12 = high. As far as I understood from other threads, the way to model the fact of nestedness also depends on how I have coded the nested factor (??). I was wondering whether I need to include the term (1|frequency/task) to specify that task is nested within frequency and whether this term is correct after all or whether it is sufficient to just include (1|task), because the task number already implies the level of frequency. Ideally, I would have liked to include task or frequency/task as a fixed effect, but the number of parameters is getting too high then. I was playing around with (1|frequency/task) and found out that I can't interpret the output either. R gives me a row with "frequency:task" and "frequency". Looks wrong, since I don't want to have frequency as a random factor... I'd be very grateful for any help. Thank you very much in advance. Eva


title: need help understanding Dirichlet (coursera's PGM class week 7 - Bayesian prediction)

text: I'm trying to work through Coursera's probabilistic graphical models class (week 7: Baeysian prediction) and a have several questions.   1. In the Dirichlet distribution, I'm having difficulty trying to understand why there's a   -1   in theta's exponent: $$P(\theta)=Dir(\alpha_1, ..., \alpha_k) = \frac{1}{Z} \cdot \prod_{j} \theta_{j}^{\alpha_{j}-1}$$   2. How do you get from here: $$P(X)=\int_{\theta}P(X|\theta)P(\theta)d\theta$$ to here: $$P(X=x^{i}|\theta) = \int_{\theta} \frac{1}{Z} \cdot \theta_{i} \prod_{j} \theta_{j}^{\alpha_{j}-1}$$   3. Also, how do you step through the integration for the following?: $$\int_{\theta} \frac{1}{Z} \cdot \theta_{i} \prod_{j} \theta_{j}^{\alpha_{j}-1} = { \alpha_{i}\over{\sum_{j} \alpha_{j}} }$$ These are the lecture notes. My questions refer to the first slide.


title: Power analysis for biased coin in R

text: I have an experiment where I'd like to determine if a coin is biased. I'm using the binomial distribution to look at coin flips. I've collected data from a sample of 30 flips and found that in 24 of the flips, the outcome was heads. I would like to calculate the power of this experiment. Here I am interested in only a one sided test:               binom.test(24, 30, p=0.5, alternative="greater")      How can I calculate the power of this analysis?


title: Transformations used in Automatic Linear Modeling procedure in SPSS

text: I am using Automatic Linear Modeling in SPSS. Does anyone here know what are the transformations option which is suggested in that procedure?


title: log-linked GLM Gamma vs log-linked GLM Gaussian vs log-transformed LM

text: From my results, it appears that GLM Gamma meets most assumptions, but is it a worthwhile improvement over the log-transformed LM? Most literature I've found deals with Poisson or Binomial GLMs. I found the article EVALUATION OF GENERALIZED LINEAR MODEL ASSUMPTIONS USING RANDOMIZATION very useful, but it lacks the actual plots used to make a decision. Hopefully someone with experience can point me in right direction. I want to model the distribution of my response variable T, whose distribution is plotted below. As you can see, it is positive skewness:   ![Valid XHTML](http://i.imgur.com/EJLgGdy.jpg). I have two categorical factors to consider: METH and CASEPART.   Note that this study is mainly exploratory, essentially serving as a pilot study before theorizing a model and performing DoE around it. I have the following models in R, with their diagnostic plots: `   LM.LOG<-lm(log10(T)~factor(METH)+factor(CASEPART),data=tdat)`   ![Valid XHTML](http://i.imgur.com/jAYBWS7.jpg)   ![Valid XHTML](http://i.imgur.com/CekkIN2.jpg) `   GLM.GAMMA<-glm(T~factor(METH)*factor(CASEPART),data=tdat,family="Gamma"(link='log')) `   ![Valid XHTML](http://i.imgur.com/4VZOa5S.jpg)   ![Valid XHTML](http://i.imgur.com/PQKcfLx.jpg) `   GLM.GAUS<-glm(T~factor(METH)*factor(CASEPART),data=tdat,family="gaussian"(link='log'))` ![Valid XHTML](http://i.imgur.com/iY9PHeh.jpg)   ![Valid XHTML](http://i.imgur.com/HV5idKy.jpg) I also attained the following P-values via Shapiro-Wilks test on residuals:   ` LM.LOG: 2.347e-11   GLM.GAMMA: 0.6288   GLM.GAUS: 0.6288   ` I calculated AIC and BIC values, but if I am correct, they don't tell me much due to different families in the GLMs/LM. Also, I noted the extreme values, but I cannot classify them as outliers as there is no clear "special cause". Thanks for any insight!


title: Possible Paradox: Calculating a confidence interval with within-experiment error

text: This is a spinoff of How to calculate the confidence interval of the mean of means? and related to When making inferences about group means, are credible Intervals sensitive to within-subject variance while confidence intervals are not? Dataset 1 here is taken from the first link above. Dataset 2 has the approximately the same experimental means but different within experiment variance. My first question is: 1) How do I calculate a confidence interval for the overall mean for each of these data sets? If I understand @Stéphane Laurent's answer in the two linked questions they should be the same.If that is true this goes strongly against all my scientific intuition and also appears to be a paradox. 2) How can it be that the confidence interval is apparently both sensitive to and not sensitive to within experiment error? ![enter image description here](http://i.stack.imgur.com/zvYXu.png) dataset 1:                 Experiment Value               1    34               1    41               1    39               2    45               2    51               2    52               3    29               3    31               3    35          structure(list(Experiment = structure(c(1L, 1L, 1L, 2L, 2L, 2L,      3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor"), Value = c(34,      41, 39, 45, 51, 52, 29, 31, 35)), .Names = c("Experiment", "Value"     ), row.names = c(NA, -9L), class = "data.frame")      dataset2:                 Experiment    Value               1 38.20744               1 37.99410               1 37.96299               2 49.27085               2 49.40519               2 49.24894               3 31.81259               3 31.73708               3 31.73834          structure(list(Experiment = structure(c(1L, 1L, 1L, 2L, 2L, 2L,      3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor"), Value = c(38.2074373061779,      37.9941025108851, 37.9629896019425, 49.2708491636015, 49.4051867974062,      49.2489418702291, 31.8125943239769, 31.7370826901692, 31.7383364604132     )), .Names = c("Experiment", "Value"), row.names = c(NA, -9L), class = "data.frame")


title: Detrending Discrete Data

text: I am trying to detrend some discrete data and I am having difficulty finding a model to describe the trend. There is a number of discrete data points and there is a linear error being introduced with time. I created some representative data in excel: Unbiased: ![Unbiased](http://i.stack.imgur.com/4fv72.png) Biased data: ![Biased](http://i.stack.imgur.com/02Qit.png) The discrete data here is 1,2,3 with a positive linear bias over time. Traditional least squares does a poor job, orthogonal least squares does no better. I tried a few other models with no success. I have been able to visually inspect the trend and select some of the ~2s, perform a linear regression, and then subtract the bias. I would like to programmatically do this across many data sets, so a generic way to do this reliability is essential. Any ideas on how to detrend the data and remove the bias? Actual data: data


title: tensor cross covariance calculation

text: I have an $N \times P$ matrix $X$ (with $\Sigma$ variance covariance matrix $P \times P$) and I am using R to calculate the tensor product $E\\{(X-\mu)(X-\mu)' \otimes (X-\mu)(X-\mu)'\\}$ (tensor product `%o%` in R). I am not sure how to evaluate this expression but I know that the result should be a $P^2 \times P^2$ matrix. Is it correct to calculate simply the tensor product using the variance covariance matrix as $\Sigma \otimes \Sigma$?


title: How to set $\nu$ parameter for one-clss SVM?

text: I am trying to perform one-class SVM for novelty detection, so I use cross- validation to find best $\nu$ parameter for it, but I found that usually I got a big parameter $\nu$, resulting in a big false positive rate. How to set and find the $\nu$ parameter for a one-class SVM?


title: Multiple imputation for missing count data in a time series from a panel study

text: I am trying to tackle a problem which deals with the imputation of missing data from a panel data study(Not sure if I am using 'panel data study' correctly - as I learned it today.) I have total death count data for years 2003 to 2009, all the months, male & female, for 8 different districts and for 4 age groups. The dataframe looks something like this:                        District  Gender Year Month    AgeGroup TotalDeaths              Northern    Male 2006    11        01-4           0              Northern    Male 2006    11       05-14           1              Northern    Male 2006    11         15+          83              Northern    Male 2006    12           0           3              Northern    Male 2006    12        01-4           0              Northern    Male 2006    12       05-14           0              Northern    Male 2006    12         15+         106              Southern  Female 2003     1           0           6              Southern  Female 2003     1        01-4           0              Southern  Female 2003     1       05-14           3              Southern  Female 2003     1         15+         136              Southern  Female 2003     2           0           6              Southern  Female 2003     2        01-4           0              Southern  Female 2003     2       05-14           1              Southern  Female 2003     2         15+         111              Southern  Female 2003     3           0           2              Southern  Female 2003     3        01-4           0              Southern  Female 2003     3       05-14           1              Southern  Female 2003     3         15+         141              Southern  Female 2003     4           0           4      For the 10 months spread over 2007 and 2008 some of the total deaths from all districts were not recorded. I am trying to estimate these missing value through a multiple imputation method. Either using Generalized Linear Models or SARIMA models. My biggest issue is the use of software and the coding. I asked a question on Stackoverflow, where I want to extract the data into smaller groups such as this:                        District  Gender Year Month    AgeGroup TotalDeaths              Northern    Male 2003     1        01-4           0              Northern    Male 2003     2        01-4           1              Northern    Male 2003     3        01-4           0              Northern    Male 2003     4        01-4           3              Northern    Male 2003     5        01-4           4              Northern    Male 2003     6        01-4           6              Northern    Male 2003     7        01-4           5              Northern    Male 2003     8        01-4           0              Northern    Male 2003     9        01-4           1              Northern    Male 2003    10        01-4           2              Northern    Male 2003    11        01-4           0              Northern    Male 2003    12        01-4           1              Northern    Male 2004     1        01-4           1              Northern    Male 2004     2        01-4           0      Going to                        Northern    Male 2006    11        01-4           0              Northern    Male 2006    12        01-4           0      But someone suggested I should rather bring my question here - perhaps ask for a direction? Currently I am unable to enter this data as a proper time- series/panel study into R. My eventual aim is to use this data and the `amelia2` package with its functions to impute for missing `TotalDeaths` for certain months in 2007 and 2008, where the data is missing. Any help, how to do this and perhaps suggestions on how to tackle this problem would be gratefully appreciated. If this helps, I am trying to follow a similar approach to what Clint Roberts did in his PhD Thesis. **EDIT:** After creating the 'time' and 'group' variable as suggested by @Matt:               > head(dat)          District Gender Year Month AgeGroup Unnatural Natural Total time                    group     1 Khayelitsha Female 2001     1        0         0       6     6    1     Khayelitsha.Female.0     2 Khayelitsha Female 2001     1     01-4         1       3     4    1  Khayelitsha.Female.01-4     3 Khayelitsha Female 2001     1    05-14         0       0     0    1 Khayelitsha.Female.05-14     4 Khayelitsha Female 2001     1     15up         8      73    81    1  Khayelitsha.Female.15up     5 Khayelitsha Female 2001     2        0         2       9    11    2     Khayelitsha.Female.0     6 Khayelitsha Female 2001     2     01-4         0       2     2    2  Khayelitsha.Female.01-4      As you notice, there's actually further detail 'Natural' and 'Unnatural'.


title: Where do the assumptions for linear regression come from?

text: I'v already known that there are several assumpations when using linear regression model. But I cannot understand why some of them exists. They are:   1. independent errors   2. normal distribution of errors   3. homoscedasticity Why cannot I just simply use least square method without these assumptions ? I want to know how it affect the $R^2$, slope, and p-value if some assumptions are invalid. For example: "if the independent error assumption is false, the p-value is less then its true value."


title: Hypothesis testing of normal distribution, known mean unknown variance

text: I've been working on review problems, and this one has me _completely_ stumped. Let $X_1 ... X_{10}$ be a random sample from a $N(3,\sigma^2)$ distribution, where $\sigma^2$ is unknown. Using the likelihood ratio test, determine a 5%-level critical region test for $H_0 : \sigma^2 = 1 $ vs. $H_1 : \sigma^2 \neq 1$ (and, trivially, $\sigma^2 >0$). It appears that in the general case, when one is testing a hypothesis about the variance, a chi-square statistic is used, which gives me something of an end-goal, but I'm not sure how to get there. The joint pdf for the 10 r.v.s should be $\large(\frac{1}{\sqrt{2\pi\sigma^2}})^{10}\cdot e^-\frac{\sum_{i=1}^{10} (X_i - 3)^2}{2\sigma^2}$ Under the null hypothesis, this yields $\large(\frac{1}{\sqrt{2\pi}})^{10}\cdot e^-\frac{\sum_{i=1}^{10} (X_i - 3)^2}{2}$, since $\sigma^2 = 1$ Under the alternative hypothesis, we have $\large(\frac{1}{\sqrt{2\pi\hat\sigma^2}})^{10}\cdot e^-\frac{\sum_{i=1}^{10} (X_i - 3)^2}{2\hat\sigma^2}$ Setting these as numerator and denominator, respectively, I get $\LARGE\frac{\exp(^-\frac{\sum_{i=1}^{10} (X_i - 3)^2}{2})}{(\frac{1}{\hat{\sigma}})^{10}\cdot \exp(^-\frac{\sum_{i=1}^{10} (X_i - 3)^2}{2\hat\sigma^2})} = \Lambda$ I believe the numerator has 0 free parameters, and the denominator has 1. In order to get the log-likelihood, I apply $ln(\Lambda)$, and we know that $\hat\sigma^2$ can be represented as $\frac{1}{10}\sum_{i=1}^{10} (X_i-3)^2$, so further simplification yields $-2Ln(\Lambda) = \sum(X_i-3)^2-10+10ln(\frac{10}{\sum(Xi-3)^2})$ According to the problem, this should be a $\chi_{10}^2$ statistic, but I don't know how to justify this (probably graphically)? Again, I greatly appreciate the help! Edit (and my proposed answer): If I instead put everything in terms of $\hat\sigma^2$, I end up with the following: $10(\hat\sigma^2-ln(\sigma^2)-1)$, and since I'm purely looking to see if this monotonic, I can simplify this to $\hat\sigma^2-ln(\hat\sigma^2)$, which a quick graph shows to be not-monotonic. This means we are going to do a two sided test under the null hypothesis. We know $\hat\sigma^2$ follows a $\chi^2_10$ distribution so we reject $H_0$ at when $n*\sigma^2< $$\chi_{.025,10}^2$ and at $n*\sigma^2>\chi_{.975,10}^2$


Suppose $X_1, \dots , X_n$ is a random sample from $f(x|\mu,\sigma)=\displaystyle\frac{1}{\sigma}\displaystyle e^{\displaystyle-\frac{x-\mu}{\sigma}}, x>\mu , \sigma>0$. How can calculate $\text{cov}(\bar{X_n},X_{(1)})$. _Note:_ $X_{(1)}$ is first order statistics.


title: Cox model vs logistic regression

text: Let's say we are given the following problem: Predict which clients are most likely to stop buying in our shop in next 3 months.   For each client we know the month when one started to buy in our shop and additionally we have many behavioral features in monthly aggregates. The 'eldest' client has been buying for fifty months; let's denote the time since a client began to buy by $t$ ($t \in [0, 50]$). It can be assumed that the number of clients is very large. If a client stops buying for three months and then comes back, then he is treated as a new customer so an event (stop buying) can occur only once. Two solutions come into my mind: **Logistic regression** \- For each client and each month (maybe except the 3 newest months), we can say whether a client stopped buying or not, so we can do rolling samples with one observation per client and month. We can use the number of months since beginning as a categorical variable to obtain some equivalent of base hazard function. **Extended Cox model** \- This problem can be also modeled using the extended Cox model. It seems that this problem is more suited to survival analysis. **Question:** What are the advantages of survival analysis in similar problems? The survival analysis was invented for some reason, so there must be some serious advantage. My knowledge in survival analysis is not very deep and I think that most potential advantages of the Cox model can also be achieved using logistic regression.   * Equivalent of stratified Cox model can be obtained using an interaction of $t$ and the stratifying variable.    * Interaction Cox model can be obtained by diving the population into several sub-populations and estimating LR for every sub-population. The only advantage I see is that Cox model is more flexible; for example, we can easily calculate the probability that a client will stop buying in 6 months.


title: Kernel Selection

text: I am not an expert in SVM and kernel, so please excuse me if I ask stupid question. Actually, first I want to know how to analyze a dataset to discover its pattern. And second, how can I select best kernel function for classifying a dataset? Thanks in advance for any help.


title: Tests for univariate outliers: have Dixon's and Grubb's methods been discredited?

text: In contrast to the many threads on this site that recommend Dixon's and Grubb's tests, the author of one answer, at this thread, contends that "Really, these have been discredited long ago" and advocates 2 other methods. I don't feel qualified to sort out these arguments, but I would like to ask whether there is any consensus among statisticians as to the merits of either position.


title: Strange: Partial Correlation changed correlations values to be negative?

text: I am trying to investigate the correlation values. A big problem happens. The method I used is in **corpcor** package in R. I have a column as missing data. I firstly calculate the correlation matrix.               cluster1corMatrix<-cor(cluster1Analytes,method="spearman")      Then I get a correlation matrix. when I did partial correlations using cor2pcor functions as follows, I got a big problem.               ## 80th column and 80th row are all NAs.     cluster1corMatrix<-cluster1corMatrix[-80,]     cluster1corMatrix<-cluster1corMatrix[,-80]     ### the above delete the missing value column     cluster1PcorMatrix<-cor2pcor(cluster1corMatrix)      Most of the correlation become negative. More importantly, some correlations, for instance, 0.9 correlation becomes -0.91 in partial correlation. The method I used is in **corpcor** package in R. Can anyone comment on that ? Is it reasonable to change correlation of 0.9 to -0.9 in partial correlation ? Many thanks for reading this post.


title: Appropriate df for a linear mixed model when looking at variance

text: Suppose I fit a linear mixed model as:               lme(Response ~ 1 , random = ~1| Location | User | Machine).      Thus machine is nested within user which is nested within location. Suppose there are 4 Locations, 8 users and 16 Machine. Suppose in total there were 96 total experiments (over all locations, users and machines) As far as I see it, I have the model $response = \mu + \alpha_i + \beta_{ij} + \gamma_{ijk} + \epsilon_{ijkl}$ where $\alpha$ corresponds to location, $\beta$ to user and $\gamma$ to machine. Thus $i \in \\{1,2,3,4\\}$, $j \in \\{1,...,8\\}$ and $k \in \\{1,...,16\\}$. Now, suppose I wanted to look at the variability of $\alpha$ and do some sort of chi squared test on it's variance. I'm assuming the df I would use would be something on the order of (4-1), since the sample size is "4". Granted, if the statistic is more complex that might change but the "n" is 4 for that. Now, if i were to do the same thing for $\beta$, I would use an n of 8. For $\gamma$ I would use an n of 16. However, what about if I was looking at the variance of $\epsilon$. Would I use an n of (96-4-8-16) = 68? Thus if my test had n-1 df I would use df = 67? Thanks!


title: model specification problem

text: I have the following model that I would like to rebuild: $Y_{i,t}=a+bx_{i,t}+cx_{i,t-1}+e_{i,t}$ I' am wondering now whether this is the same as the model above: $Y_{i,t}=a+ d\Delta x_{i,t}+e_{i,t}$, where $\Delta x_{i,t}=x_{i,t} - x_{i,t-1}$? I believe this is true but i just need to be 100% sure. Thank you


title: Combining two interval measures to correlate to a third - How?

text: I want to compare the scores received on two separate measures together to a third measure. They each have raw scores available and I was wondering if I found the z-scores for the first two interval measures and combine/average could I then compare thre results to a third interval measure for correlation?


title: Linear Factor Model vs. Linear Regression Model

text: I've been reading some literature that discusses 'linear factor models' which appear to describe the general equation often used in OLS regression. When people refer to a 'linear regression model' are they essentially just referring to a linear factor model? Where does the term linear factor model fit in in statistics and why, if at all, is it necessary to make the distinction?


title: 3 Conceptual Questions about a DBN (Deep Belief Network)

text: Paper is from here: https://www.cs.toronto.edu/~hinton/absps/fastnc.pdf (1) Why is there a RBM sitting on top? (2) Why is a RBM considered equivalent to an infinite sigmoid network? (3) Is a chain of stacked RBMs a DBM (Deep Boltzmann Machine) or a DBN?


title: Difference in output from mvtnorm and Mathematica/Java

text: Let $a \sim\mathcal{N}(6.532056,0.06532056)$,$b \sim\mathcal{N}(8.390961,0.08390961)$ and $c \sim\mathcal{N}(8.736566,0.08736566)$. We use $\mathcal{N}(\mu,\sigma^2)$ notation unless specified otherwise. We construct two normal variables $x = a-b$ and $y = a-c$. So, $x \sim\mathcal{N}(-1.858905,0.14923017)$ and $y \sim\mathcal{N}(-2.20451,0.15268622)$ Correlation between $x$ and $y$ i.e. $\textrm{cor}(x,y)= \frac{\textrm{var}(a)}{\sigma_x \times \sigma_y}$ by using the basic properties of covariance. Solving, $\textrm{cor}(x,y)=0.4327346392418512 \approx 0.433$ which can be written in matrix (let it be called $ mat $)for as \begin{bmatrix} 1.0 & 0.433 \\\\[0.3em] 0.433 & 1.0 \\\\[0.3em] \end{bmatrix} Now I want to find the $\Pr[-\infty<x<0 \textrm{ and } -\infty<y<0]$. I used the `mvtnorm` package's `pmvnorm` method by invoking               pr<-pmvnorm(mean=c(-1.858905,-2.20451), corr=mat, lower=rep(-Inf, 2), upper=rep(0,2))      The result was 0.9575448. The same when I want to compute in `Mathematica/Java`(an implementation of the Genz algorithm found on internet) I am getting the result as 0.9999992445813132. I am including the Mathematica code below. Please note here, the second argument in this case, is standard deviation.               px = NormalDistribution[-1.858905000000001, 0.38630320992712447]          py = NormalDistribution[-2.20451, 0.3907508413298685]          pxy = ProductDistribution[px, py]          Probability[-Infinity < x < 0 && -Infinity < y < 0, {x,         y} \[Distributed] pxy]      What is the mistake I am doing? For my application domain (visualization) this difference is turning out to be very costly,especially because I run it over more than 700,000 data points.


title: logistic regression using probabilities of class labels

text: My goal is to train a logistic classifier. My samples in my dataset have some label noise but for each label I can give a probability how correct this label is. What is the best way to incorporate these probabilities while fitting my model? I'm using the glmnet package in R for my experiments. According to the help you can give a 2-column vector as response variable. So is it correct to give the probabilities of being class 1 in the first column, probabilities of being class 2 in the second column. And set the weights to 1 for every sample. Is this correct?


title: Multiplicative error and additive error for generalized linear model

text: If the following generalized linear model was used, how should I interpret the error term?   link function: natural log   distribution: Gamma distribution   i.e., $\ln E(Y)=X\beta$ and $E(Y)=\exp(X\beta)$   It seems that the error term should be additive:   $Y=\exp(X\beta)+\epsilon$   However, after some algebra, it becomes:   $Y=\exp(X\beta)\psi$ by letting $\psi_i=\frac{\epsilon_i}{\exp(x_i'\beta)}+1$   Thus, $\epsilon$ is additive but $\psi$ is multiplicative.      I am confused with this article:   "Multiplicative errors: Log-normal or gamma?" by Firth, 1988   Firth established the multiplicative errors model for gamma distributed response variable.   However, after the computation by R or SAS, how do I interpret the error term?   Does R glm() or SAS proc genmod provide us the multiplicative error model (as in Firth, 1988)?      Moreover, why is the following model not available in R or SAS?   distribution: lognormal distribution   link function: identity   It seems (in many literatures) that we have to log-transform the response, and apply the following model (or using OLS):   distribution: normal distribution   link function: identity   However, in comparison to the gamma model (log link), we must construct the following model:   distribution: lognormal distribution   link function: natural log   But it could not be estimated in R glm or SAS proc genmod.   After some research on the topic "natural exponential family",   http://en.wikipedia.org/wiki/Natural_exponential_family   I guess that the algorithm in R glm() or SAS proc genmod cannot solve models in which the distribution is not in the natural exponential family. (if the Fisher-scoring is used, not quasi-likelihood)   It seems that the "quasi-likelihood" must be applied in order to solve the models for lognormal models without log-transform.   However, these are all my guess. I need some literatures (journal articles, books, etc.) to be cited in order to convice others.   Or a more rigorous proof may be helpful.


title: Confusion related to matlab implementation of liblinear

text: I have this confusion related to the matlab implementation of liblinear. I am trying to use L1 regularized logistic regression. I have trained the model with a bias. Now when I do the prediction calling               [predicted_label, accuracy, decision_values/prob_estimates] = ...     predict(testing_label_vector, testing_instance_matrix, model, 'liblinear_options','col']      Do I need to add an extra feature with ones for the bias term. I mean when I learn the model if I have 10 features, it will have 11 weights including the one for the bias Now when I provide my dataset for testing, do I need to add an extra column in the `testing_instance_matrix`


title: Alternatives to proportional odds mixed effects models

text: I was doing some reading about categorical data analysis (e.g. Agresti (2010): Analysis of Ordinal Categorical Data). Additionally to analysing data with an ordinal response with the help of proportional odds cummulative mixed effects models it seems some people suggest e.g. to use the `glmer()` function from `lme4` in `R` to fit a whole range of models. My data looks as follows:   1. ordinal response: a five point scale with $1<2<3<4<5$ ($1$ being worst, $3$ being neutral, $5$ being best)   2. one single categorical predictor "type" with 5 levels   3. subjects as random effects term which I would specify in `lme4` syntax: `(type | subjects)` with intercept and slope for every subject. Every subject saw $320$ items ($40$ for level $1$, $40$ for level $2$, $40$ for level $3$, $160$ for level $4$) and rated every single one of them. The dataset seems quite nice and as it only contains a single predictor is very good to get acquainted with categorical data analysis. However, I was wondering if there are alternative ways of analysing such data apart from treating the response as a continuous variable and using `lmer()` or `aov()`. And if there are any recommendation or advice people can give here. Here is an alternative I read about . One could split a dataset with an ordinal response into subsets that each only contain two levels of the ordinal response. There are some questions I have about this. How would you split that dataset into subsets given a $5$ point ordinal response? Based on the ordering of my response the most reasonable choice to me seems: subset $A$ with levels $1$ and $2$, subset $B$ with levels $2$ and $3$, subset $C$ with levels $3$ and $4$, subset $D$ with levels $4$ and $5$. I tried this with the `wine` data set from the `ordinal` package. But this will ungraciously fail because of convergence failures, model unidentifiability etc. (I'm using `lme4 1.1.8`). Here is the wine example:               data("wine", package = "ordinal") # the wine dataset     library(lme4)     sub <- wine$rating %in% grep("[1-2]", wine$rating, value = TRUE)     sub12 <- subset(wine, subset = sub)     sub12$rating <- factor(sub12$rating, ordered = FALSE)     glmer12 <- glmer(rating ~ temp + contact + (temp | judge),                      data = sub12,                      family = binomial)          # Substituting (1 | judge) for (temp | judge) won't make it better          # Will give no warning     glm12 <- glm(rating ~ temp + contact,                  data = sub12,                  family = binomial)          sub <- wine$rating %in% grep("[2-3]", wine$rating, value = TRUE)     sub23 <- subset(wine, subset = sub)     sub23$rating <- factor(sub23$rating, ordered = FALSE)          glmer23 <- glmer(rating ~ temp + contact + (temp | judge),                      data = sub23,                      family = binomial)          sub <- wine$rating %in% grep("[3-4]", wine$rating, value = TRUE)     sub34 <- subset(wine, subset = sub)     sub34$rating <- factor(sub34$rating, ordered = FALSE)     glmer34 <- glmer(rating ~ temp + contact + (temp | judge),                      data = sub34,                      family = binomial)      Of course the wine dataset is rather tiny (and we have full- or quasi- separation) and mine won't but how likely am I to run in similar problem if I analyse ordinal data this way?


title: lme() with several within and between (categorical and continuous) subject factors

text: I am currently trying to analyse data from an experiment of mine and I have done some searching for instructions on the usage of the lme() function for R, since I am looking to analyse my data with a linear mixed effects approach. However,   1. In case lme() does work out, I am not sure how to fit a model to my data.   2. In case there are better options than lme(), I would be glad about any help regarding the proper use of these functions. Subjects in my experiment were randomly assigned to one of two contexts (between subjects factor B) and within this context all of the subjects completed a task which was constructed as a combination of two within subjects factors (W1 and W2). Additionally, subjects have filled a personality questionnaire, thus each subjects has a value for this, which I need to include in my analysis (Pers). The outcome variable is a continuous one (reaction time, RT). The data look like this:               B   Pers    subject W1  W2  RT     A   23      1       m   x   187     A   23      1       m   y   333     A   23      1       n   x   112     A   23      1       n   y   110     A   5       2       m   x   313     ... ...     ...     ... ... ...     B   11      19      m   x   911     B   11      19      m   y   248     B   11      19      n   x   411     B   11      19      n   y   212     B   -9      20      m   x   666     ... ...     ...     ... ... ...      I am interested in an analysis giving me information about the interactions of all four factor (B, Pers, W1 and W2), thus four-way interaction. When including only B, it seems to me that I understand what to do when using lme(). My code looks like this (nb: I constructed contrasts between the different conditions before the actual lme() code using the contrasts() function):               > baselinemodel = lme(RT ~ 1, random = ~1|subject/W1/W2, data = df, method = "ML")     > W1model = update(baseline, .~. +W1)     > W2model = update(W1model, .~. +W2)     > Bmodel = update(W2model, .~. +B)     > W1W2mod = update(Bmodel, .~. +W1:W2)     > W1Bmod = update(W1W2mod, .~. +W1:B)     > W2Bmod = update(W1Bmod, .~. +W2:B)     > fullmodel = update(W2bmodel, .~. +W1:W2:B)      I.e. I set up a basline model including only the intercept and specifying W1 and W2 as nested within subjects (hence defining them as within subjects factors) and then in a sequential manner I include more predictors and their interactions in order to finally be able to compare these models using the anova() function:               > anova(baseline, W1model, W2model, Bmodel, W1W2mod, W1Bmod, W2Bmod, fullmodel)      Naturally I also take a look at my predefined contrasts using               > summary(fullmodel)      I find this approach very intuitive, as long as I do not try to include the continuous predictor Pers in addition to B. So this is where I am stuck. I would be deeply grateful if anyone could   1. explain to me how to include Pers in addition to B, W1 and W2 (i.e. how to include a continuous predictor varying between subjects, and its interactions with the other factors into my analysis)   2. and in case lme() is not the optimal function for my objectives how to do a comparable analysis with another function (thus broadening my still quite narrow scope of R) I have, naturally done some searching before starting to post this question, but unfortunately others, who seemed to have related issues weren't too successful in getting answers (Analysis of a mixed design with categorical and continuous variables) Thank you in advance


title: Using k-means with other metrics

text: So I realize this has been asked before: e.g. What are the use cases related to cluster analysis of different distance metrics? but I've found the answers somewhat contradictory to what is suggested should be possible in the literature. Recently I have read two papers that have mention using the kmeans algorithm with other metrics, for example edit distance between strings and the "Earth Mover Distance" between distributions. Given that these papers mention using kmeans with other metrics without specifying _how_ , particularly when it comes to computing the mean of set of points, suggests to me that maybe there is some "standard" method to dealing with this that I'm just not picking up on. Take for example this paper, which gives a faster implementation of the k-means algorithm. Quoting from paragraph 4 in the intro the author says his algorithm "can be used with any black box distance metric", and in the next paragraph he mentions edit distance as a specific example. His algorithm however still computes the mean of a set of points and doesn't mention how this might affect results with other metrics (I'm especially perplexed as to how mean would work with edit distance). This other paper describes using k-means to cluster poker hands for a texas hold-em abstraction. If you jump to page 2 bottom of lefthand column the author's write "and then k-means is used to compute an abstraction with the desired number of clusters using the Earth Mover Distance between each pair of histograms as the distance metric". I'm not really looking for someone to explain these papers to me, but am I missing some standard method for using k-means with other metrics? Standard averaging with the earth mover distance seems like it could work heuristically, but edit distance seems to not fit the mold at all. I appreciate any insight someone could give. **(edit)** : I went ahead and tried k-means on distribution histograms using the earth mover distance (similar to what is in the poker paper) and it seemed to have worked fine, the clusters it output looked pretty good for my use case. For averaging I just treated the histograms as vectors and averaged in the normal way. The one thing that I noticed is the sum over all points of the distances to the means did not always decrease in a monotone manner. In practice though, it would settle on a local min within 10 iterations despite monotone issues. I'm going to assume that this is what they did in the second paper, the only question that remains then is, how the heck would you average when using something like edit distance?


title: Why normalise the standard deviation of neural network input?

text: I understand why the mean of the input to a neural network is normalized, to avoid the numerical problems with very large and very small numbers. Also, it's nice if the bias node is around the same magnitude as the other input data. But why is the standard deviation usually normalized? Of course higher standard deviation means that the variations have larger effect, but won't the weights adapt to that anyway? Thanks!


title: What prior distributions are used in mcmcsamp() from lme4?

text: The `mcmcsamp()` function generates simulations from the posterior distributions of a Bayesian mixed model fitted with the `lmer()` function. What are the prior distributions? Does there exist a book or an article providing a mathematical study of the Bayesian mixed model with these priors?


title: Goodness of fit between empirical distribution and theoretical distribution

text: Suppose one wants to assess the goodness of fit between a theoretical distribution and an empirical distribution in graphical manner. The normal Q-Q plot is good if the theoretical distribution is normal. But what if the theoretical distribution is not normal? Can you use a double logarithmic graphical method?


title: Effect of each parameter on a Monte Carlo Simulation

text: I was wondering what is the best way to determine the effect of each random parameter on the result obtained from a _Monte Carlo Simulation_. I realise I have asked a similar question here, but this time I am making this question in a more general manner as I didn't get any answers previously. **Background** : I'm essentially running a Monte Carlo Simulation where I sample randomly from 65 Normal Distributions and 2 Uniform Distributions. Then, i then input these values to a linear equation and run this for 500.000 iterations. I just wanted to determine the best way to establish the relative effect of each of these 67 (in total) random variables on the final outcome from my Monte Carlo Simulation.


title: Bootstrapping in R

text: Is there anyway to perform a BCa, a tilted, and an ABC bootstrap in R without using the built in packages, namely 'boot'? Thanks


title: How to identify outliers from a Gumbel distribution with known parameters?

text: I have a model for data in my experiment that states that the data has a Gumbel distribution with known location and scale. I am then looking at observations with very high scores that I suspect to be outliers. Is there a Bayesian approach (or otherwise) to get the probability that this observation is from the underlying Gumbel distribution? As an end result I would like to decide whether or not I accept this observation as coming from the underlying Gumbel distribution. ### Edit I am in fact working with biological data and have demonstrated both empirically and through a set of assumptions that the distribution is in fact of Gumbel type with parameters that can be estimated from the data itself and/or from a theoretical perspective given some assumptions. Now, while most data behaves in-line with this distribution, I would like to find the probability that observations with very large values might _not_ be from the same Gumbel distribution but due to some other phenomenon that is not described in this distribution that I want to detect. In this case, updating the data using the new observation would not be useful (especially considering that this observation might come from another random variable with a different distribution). I am able to calculate the one-side probability of how extreme the observation is by using. $$1-F(x)$$ From my (very limited) understanding of statistics and probability, I fear that quoting this value might be erroneous. My main concern is that this is the "probability of observing more extreme observations" and _NOT_ the "probability that this observation is not from the distribution". Any thoughts on this?


title: Marginal prior distribution with restricted parameters

text: I'm analysing one paper on bayesian inference for network reliability and got stuck at trying to validate some (quite simple at first sight) formulas. Suppose the probabilities of failure has truncated beta distribution as a prior:   $$ \pi _{1}\left ( p_{i}|\gamma ,\alpha \right )=\frac{\alpha \left ( 1-p_{i} \right )^{\alpha -1}}{\left ( 1-\gamma \right )^{\alpha }}, i=\overline{1,n}, \gamma <p_{i}<1$$ where $\gamma$ is a lower threshold on $p_{i}$ and $\alpha$ is just a fixed parameter (i.e. has no prior density). The second stage prior is for parameter $\gamma$ and is a beta distribution with parameters $ \alpha +1 $ and $1$. That is,   $$\pi _{2}\left ( \gamma |\alpha \right )=(\alpha +1)(1-\gamma)^{\alpha},\, 0<\gamma <1, \, 0<\alpha<1$$ I want to find prior of $p_{i}$ conditional on $\alpha$ alone, i.e. $\pi\left( p_{i}|\alpha\right)$. I apply the decomposition $$\pi \left ( p_{i}, \gamma|\alpha \right )=\pi_{1}\left ( p_{i}|\gamma,\alpha \right )\cdot \pi_{2}\left ( \gamma|\alpha \right )$$ And marginalization over $\gamma$ should give what I'm looking for: $$\pi\left( p_{i}|\alpha\right)=\int_{0}^{1}\alpha\left ( \alpha +1 \right )\left ( 1-p_{i} \right )^{\alpha -1}d\gamma=\alpha\left ( \alpha +1 \right )\left ( 1-p_{i} \right )^{\alpha -1}$$ But this does not integrate to one and the paper gives another answer (which does integrate to one): $$\pi\left( p_{i}|\alpha\right)= \alpha\left ( \alpha +1 \right )p_{i}\left ( 1-p_{i} \right )^{\alpha -1}$$ And I can't figure out what I'm missing. Maybe someone could give me a clue what I'm missing?


title: How to justify a distribution fit is good

text: I've got 5 datasets to which I would like to fit a model distribution. I'd like to use the same distribution for each dataset but with different parameters. So I use MLE to compute the best parameters for a whole bunch of distributions and I take the one that "best" fits all 5 datasets. Here is where I am stuck a bit. In my first attempt, for each of the datasets and each model distribution I computed the $L_2$ norm of the error between the original dataset and the predicted value from the model. This gave me an error vector of 5 values which I took the magnitude of and called the "best" fit the one with the lowest magnitude. Doing this, I get the distribution that (to my eye) looks like the best fit (I plotted them all just to compare). This is exciting, but I don't know if it's rigorous enough. Then I decided to try another method and performed the Kolmogorov-Smirnov with my original dataset and the model distribution from the MLE. I stored the $D$ values in my error vector and took the magnitude of that vector with the minimum magnitude being the best. And it gave the same distribution as the $L_2$ approach! I was excited until I realized that the KS test shouldn't be used to test a model distribution against the data from which it was generated. So then I did some more reading and for each of the 5 original datasets used in the MLE, I generated a new dataset using boostrap resampling with replacement. I used these resampled datasets to compare against the model distributions and stored the $D$ values in my error vector and again took the magnitude, with the best fit being the lowest magnitude. This time it gave a different distribution as best. One that does look pretty good (by eye) but doesn't look nearly as good as the original (again, by eye). My experience in this kind of work is very limited but now I don't know what to trust. Both distributions are actually good fits. The one chosen by the $L_2$ norm is the one I think is best, but that's hardly rigorous. It's concerning that it's also the one chosen by the KS test with the original data because the test is flawed under those conditions. It may be possible that the $L_2$ method and the KS test with resampling are both correct and it's essentially the same as getting different answers based on the choice of norm (in this case, $L_2$ norm and probably the $L_\infty$ norm). Is something I did flawed or is there a way to determine which approach to chose as better than the other?


title: how to find a linear combination of predictors maximizing correlation between its score and dependent variable in R

text: Please correct me if I am wrong as I am not good at R. I think I can find a linear combination maximizing correlation between predictors and dependent variables by running partial least squares analysis with standardized X without intercept and standardized Y matrix. But how can I do that with pls or any other package in R? Does it do that with default settings?


title: Interpretation of the output from R's penalized linear discriminant analysis, using the penalizedLDA package

text: I am trying to use a penalized linear discriminant analysis for feature selection.The `R` package used was `penalizedLDA`. From my understanding, I would need to know the number of vectors, k, used and the "loadings" from the vectors to decide the top 10 (lets say) features. I gather I would have to select the absolute value of the 10 highest loadings. Can anyone help explain how one should "properly" interprete the output from running a penalized LDA using the `penalizedLDA` package. The structure of my output data:               > str(out)     List of 10      $ discrim: num [1:401, 1:35] -0.00842 -0.01026 -0.01266 -0.01543 -0.01879 ...          $ xproj  : num [1:6937, 1:35] 40 52.4 52.8 38.1 53.6 ...      $ K      : int 35          $ crits  :List of 35       ..$ : num [1:4] 1165 1165 1165 1165           ..$ : num [1:4] 294 294 294 294       ..$ : num [1:5] 25.9 25.9 25.9 25.9 25.9           ..$ : num [1:4] 19.7 19.7 19.7 19.7       ..$ : num [1:4] 15.9 15.9 15.9 15.9           ..$ : num [1:4] 3.85 3.85 3.85 3.85       ..$ : num [1:4] 1.88 1.89 1.89 1.89           ..$ : num [1:4] 1.35 1.35 1.35 1.35       ..$ : num [1:6] 0.801 0.801 0.801 0.801 0.801 ...           ..$ : num [1:4] 0.793 0.793 0.793 0.793       ..$ : num [1:4] 0.508 0.508 0.508 0.508           ..$ : num [1:4] 0.378 0.378 0.378 0.378       ..$ : num [1:4] 0.28 0.28 0.28 0.28           ..$ : num [1:4] 0.164 0.164 0.164 0.164       ..$ : num [1:4] 0.147 0.147 0.147 0.147           ..$ : num [1:4] 0.0943 0.0943 0.0943 0.0943       ..$ : num [1:6] 0.0906 0.0906 0.0906 0.0906 0.0906 ...           ..$ : num [1:4] 0.0774 0.0774 0.0774 0.0774       ..$ : num [1:4] 0.0528 0.0528 0.0528 0.0528           ..$ : num [1:5] 0.0442 0.0442 0.0442 0.0442 0.0442       ..$ : num [1:4] 0.037 0.037 0.037 0.037           ..$ : num [1:5] 0.0317 0.0317 0.0317 0.0317 0.0317       ..$ : num [1:4] 0.0223 0.0223 0.0223 0.0223           ..$ : num [1:4] 0.0172 0.0172 0.0172 0.0172       ..$ : num [1:4] 0.0126 0.0126 0.0126 0.0126           ..$ : num [1:4] 0.0104 0.0104 0.0104 0.0104       ..$ : num [1:4] 0.00787 0.00787 0.00787 0.00787           ..$ : num [1:4] 0.00576 0.00576 0.00576 0.00576       ..$ : num [1:4] 0.00454 0.00454 0.00454 0.00454           ..$ : num [1:4] 0.00436 0.00436 0.00436 0.00436       ..$ : num [1:4] 0.00328 0.00328 0.00328 0.00328           ..$ : num [1:4] 0.00281 0.00281 0.00281 0.00281       ..$ : num [1:4] 0.00225 0.00225 0.00225 0.00225           ..$ : num [1:4] 0.00186 0.00186 0.00186 0.00186       ..$ : num [1:4] 0.00142 0.00142 0.00142 0.00142          $ type   : chr "standard"      $ lambda : num 0.001          $ lambda2: NULL      $ wcsd.x : Named num [1:401] 0.368 0.327 0.318 0.367 0.478 ...           ..- attr(*, "names")= chr [1:401] "Spec400" "Spec405" "Spec410" "Spec415" ...          $ x      : num [1:6937, 1:401] 0.13038 -0.08583 0.00923 0.16059 0.02393 ...       ..- attr(*, "dimnames")=List of 2       .. ..$ : NULL           .. ..$ : chr [1:401] "Spec400" "Spec405" "Spec410" "Spec415" ...      $ y      : num [1:6937, 1] 1 1 1 1 1 1 1 1 1 1 ...      - attr(*, "class")= chr "penlda"


title: Generate random multivariate values from empirical data

text: I'm working on a Monte Carlo function for valuing several assets with partially correlated returns. Currently, I just generate a covariance matrix and feed to the the `rmvnorm()` function in R. (Generates correlated random values.) However, looking at the distributions of returns of an asset, it is not normally distributed. This is really a two part question:   1) How can I estimate some kind of PDF or CDF when all I have is some real- world data without a known distribution? 2) How can I generate correlated values like rmvnorm, but for this unknown (and non-normal) distribution? Thanks!


title: What if Lasso selects transformed terms but not untransformed terms

text: Suppose I have standard normal features $X_i \in \\{X_i : i \in \\{1,...,1000\\}\\}$. I extend this set of predictors with transformations as follows: $\\{X_i,X_i^2,X_iI(X_i > 0) : i \in \\{1,...,1000\\}\\}$. What happens if Lasso would pick $X_i^2$ or $X_iI(X_i > 0)$ but not $X_i$ itself. What do I do? Is this even a problem?


title: What model should one use for this short time series?

text: Below I have quarterly total sales on the left (dependent variable), and a sample of the sales on the right. The two variables share a correlation of 98.7%. What model should I use to predict X? For that model, should I include a constant? Seasonal adjustments? Remove outliers? The most important criteria is minimizing out of sample prediction error.               Q3'10   40.19   0.2386     Q4'10   39.36   0.2000     Q1'11   51.25   0.2173     Q2'11   54.99   0.2630     Q3'11   50.38   0.2242     Q4'11   50.77   0.2623     Q1'12   67.39   0.3548     Q2'12   77.14   0.3716     Q3'12   72.54   0.3451     Q4'12   80.21   0.3816     Q1'13   94.57   0.4661     Q2'13   102.13  0.4919     Q3'13E  X       0.4424


title: Incremental SVD in Collaborative Filtering

text: In the so-called incremental SVD used for collaborative filtering: http://www.machinelearning.org/proceedings/icml2007/papers/407.pdf http://www2.research.att.com/~volinsky/papers/ieeecomputer.pdf http://www.quuxlabs.com/blog/2010/09/matrix-factorization-a-simple-tutorial- and-implementation-in-python/ The user x item matrix R is factored as QP using gradient descent. In the classical SVD there is the diagonal matrix S which holds the singular values. What happens(ed) to that matrix in this formulation? Is it just omitted and they still call it SVD or is it implicitly part of Q and/or P?


title: Bayesian and frequentist interpretations vs approaches

text: I have been reading about the frequentist vs bayesian issue (this article has helped a lot, specially with the example; also this one), and I haven't come to terms with it. At the moment it seems like there are the frequentist and bayesian interpretation of probability; and, separately, the frequentist and bayesian approach to problems. The former is about the belief vs frequency issue (illustrated in the second article). The later is illustrated in the first article. Both put together seem to me like this:   * The frequentist interpretation of the frequentist approach ensures to be right a% of the time for large number of trials assuming only the likelihood distribution, no matter which parameter we get, as long we as assume that we'll get a good range of data.   * The frequentist interpretation of the bayesian approach ensures to be right a% of the time for large number of trials assuming the likelihood distribution and the prior, no matter which data we get, as long as we assume that we'll get a good range of parameters.   * The bayesian interpretation of the frequentist approach says that we are right with a probability of a% assuming only the likelihood distribution, no matter which parameter we get, as long as we assume fairness in the randomness of the data.   * The bayesian interpretation of the bayesian approach says that we are right with a probability of a% assuming the likelihood distribution and the prior, no matter which data we get, as long as we assume fairness in the randomness of the parameters. This is the only consistent view that I have been able to form from what I've read. However, I I still think I maybe missing something (as I actually haven't found this view like this anywhere else, it's my own conclusion), so taking the null hypothesis that I'm wrong, where's my mistake?


title: An instance of the Cauchy–Schwarz inequality

text: In the proof of Theorem 6.5 from the book by Devroye et al., how is the last inequality derived? $$ \begin{aligned} \mathbb{E}\left\\{|\eta(X)-1/2|\mathbf{I}_{\\{g(X)\ne g^*(X)\\}}\right\\} &\leq \mathbb{E}\left\\{\mathbf{I}_{\\{\eta(X)\ne1/2\\}}|\eta(X)-\tilde\eta(X)|\mathbf{I}_{\\{g(X)\ne g^*(X)\\}}\right\\}\\\ &= \mathbb{E}\left\\{|\eta(X)-\tilde\eta(X)|\mathbf{I}_{\\{g(X)\ne g^*(X)\\}}\mathbf{I}_{\\{|\eta(X)-1/2|\leq\epsilon\\}}\mathbf{I}_{\\{\eta(X)\ne1/2\\}}\right\\}\\\ &\+ \mathbb{E}\left\\{|\eta(X)-\tilde\eta(X)|\mathbf{I}_{\\{g(X)\ne g^*(X)\\}}\mathbf{I}_{\\{|\eta(X)-1/2|>\epsilon\\}}\right\\}\\\ &\leq \sqrt{\mathbb{E}\left\\{(\tilde\eta(X) - \eta(X))^2\right\\}}\\\ &\times \Bigg(\sqrt{\mathbb{P}\left\\{|\eta(X)-1/2|\leq\epsilon,\eta(X)\ne1/2\right\\}}\\\ &+\sqrt{\mathbb{P}\left\\{g(X)\ne g^*(X),|\eta(X)-1/2|>\epsilon\right\\}}\Bigg) \end{aligned} $$ Note that $\eta(x) = \mathbb{E}\\{Y|X=x\\}$ is the regression function, $\tilde\eta(x)$ is an approximation of $\eta(x)$, $g^*(x)$ is the Bayes classifier $$ g^*(x) = \begin{cases} 0 & \text{if } \eta(x)\leq\dfrac{1}{2} \\\ 1 & \text{otherwise} \end{cases}$$ and finally, $g(x)$ is defined like $g^*(x)$ with $\tilde\eta(x)$ replacing $\eta(x)$. $\epsilon>0$ is fixed. $\mathbf{I}_A$ is the indicator function of the set $A$.


title: Estimating an underlying pdf from binomial trials

text: I'm afraid I'm not an expert in statistics, but I have a particular problem I'm interested in solving. I'm pretty sure this area already has a lot of literature, but I'm having difficulty finding something directly applicable to what I'm doing, so it would be great if somebody could nudge me in the right direction. I'll give an example to illustrate the problem: Say I have a machine which produces biased coins. It has some underlying continuous probability density function which it is using to pick a number 0 < p < 1, then creates a coin which will come up heads with probability p. My task is to estimate the function which the machine uses to generate these coins. I'm allowed to flip the coins, and I have access to a very large number of coins, but each one after a certain, random, amount of flips is taken away from me. The number of flips is not necessarily large. How would I go about doing this? My initial thought was to just sum each resulting binomial distribution and divide by the total number of coins tested. But I'm pretty sure this doesn't give a good result. I'm faintly aware of kernal density estimation, but I don't have enough expertise to know if/how it can be used for this kind of task, or what I should know in terms of tailoring it to this task.


title: Trying out a $t$-test using R

text: I would like to: i) check if the odds ratios (OR) for each factor level are significantly different. **Edited** : This is the model I used:               NewMod2C <- glm(formula = surv ~ as.factor(season) + as.factor(pgrp) + as.factor(sline)                                                        + as.factor(tb5)                                                         + as.factor(gest3) + as.factor(int3)                                                        + as.factor(agit)                                                         + as.factor(teat2),                                                        family = binomial(link = "logit"),                                                        data = lesna)     summary(NewMod2C)          confint(NewMod2C)          data.frame(exp(NewMod2C$coefficients))                             exp.NewMod2C.coefficients.     (Intercept)                         6.5599000     as.factor(season)2                  0.7203071     as.factor(pgrp)2                    1.3495993     as.factor(pgrp)3                    1.0388424     as.factor(pgrp)4                    1.1013092     as.factor(sline)2                   1.2668896     as.factor(tb5)2                     0.8663428     as.factor(tb5)3                     0.7307868     as.factor(tb5)4                     0.6555165     as.factor(tb5)5                     0.5770979     as.factor(gest3)2                   1.3520849     as.factor(gest3)3                   2.6598991     as.factor(int3)2                    0.8146875     as.factor(int3)3                    0.2941530     as.factor(agit)2                    0.6305899     as.factor(agit)3                    0.7208870     as.factor(teat2)2                   1.3266645     as.factor(teat2)3                   1.5567743      Are there any ways of determining if, for example, the levels of `pgrp` are significantly different?                   as.factor(pgrp)2                    1.3495993 a         as.factor(pgrp)3                    1.0388424 b         as.factor(pgrp)4                    1.1013092 c      Thanks.


title: Formulating a mixed model with nested fixed and random effects

text: I'm analysing the results of an experiment using a mixed (fixed and random- effects) ANOVA with nested effects and I just need someone to confirm that I have the right model. The effects in the design are:   * Product--fixed, 5 levels--this is the main effect of interest   * Batch--random, 3 levels--nested within product   * Day--random,6 levels--this is a blocking factor in the design. On each day, multiple samples from each product/batch are tested in a random order.  So we have a complete, balanced experiment within each day. Here is the ANOVA model that I think I need:               Day      + Product     + Batch(Product)     + Day x Batch(Product)      Am I right in thinking that the appropriate error term for Product is the Batch-within-Product Mean Square? I have used type III sums of squares. I'm only really interested in testing for difference between products, so the model in the random factors is mostly only important to ensure that I get the correct residual for this comparison. However, tests for effects due to Day and Batch-within-Product are of secondary interest.


title: Edge effects in K-function

text: I am trying to write some codes for distance-based evaluation of point process such as $G$, $F$ and $K$ functions. In the implementation of $K$ there is a case to consider **edge effect** to have an unbiased $K$ measure. A simple (!) but perfect (also common) solution is to consider the intersection area of a circle drawn on each point ($center(x,y)$) with a specified radius ($r$) with the region of study. Although the concept ($K$) is purely statistical the intersection of a circle with region (in simple case a rectangle but normally a polygon) is pure geometrical problem. Any idea on:   * **how to find intersection area of a circle with arbitrary shape polygon**?   * **how to find weights for $K$ more efficiently**?


title: Several covariates in ANCOVA? Or difference scores?

text: I am running a quasi-experimental study where I compare changes in brand attitudes across time for four brands, all within subjects. Between Time 1 and time 2, participants see some information and i want to see how it impacts their rating of each brand. I want to use ANCOVA to control for Time 1 brand attitudes and then look for changes in attitudes at Time 2. Should I just put all four Time 1 brand attitudes in the ANCOVA as covariates? This seems odd since i feel like they should be matched to their respective Time 2 attitudes. DV is brand attitudes. Or should I be using difference scores? Any help would be appreciated. Thanks. Alex.


title: Maximizing incomplete likelihood

text: Given the conditional distribution $p(x|y)$ and the prior of the hidden variables $p(y|\theta)$ with unknown hyper-parameter $\theta$. Now we have observed i.i.d. samples of $x$. Besides the Bayes and empirical Bayes approach, is it possible to estimate both the hidden variables $y$ and the hyper-parameter $\theta$ directly by maximizing the incomplete likelihood $p(x,y|\theta)$? Form the Bayesian perspective, it is equivalent to finding the joint mode of $p(y,\theta|x)$ given an uninformative flat prior of $\theta$.


title: Comparison of cell counts with a right censoring

text: I have cell counts related to the action of different microorganisms and I want to compare their distribution. It's supposed they follow a normal distribution after a log transformation, but I can't register all the data due to the limit of measuring instrument, so I obtain a very similar normal curve but it's cut in the point in which I can't register the data, and if I run a K-S test, it doesn't support the normality hypothesis (but n is large). My questions: I find this can be an usual problem, so how would you demonstrate normality in these cases? or is it impossible? I think a Kruskal-Wallis test could works in my case, but I don't know if I have to take into account something else. Someone has suggested me to use an ANOVA due to my n is very high (n>2000) in each microorganism, despite the lack of normality, is this right? After reading comments I think I have a right censoring. I have seen in survival analysis is used a Kaplan-Maier estimator, but this case is different. Would it be better just remove censored data and run a non- parametric test?


title: $Cov(\hat{\epsilon})$ in a linear regression model

text: ![enter image description here](http://i.stack.imgur.com/VIoz3.jpg) I proved the result myself too, but my proof was a lot lengthier. Apparently however, it can be done in just these two lines. However, I don't really get the first step here. Could anyone please explain me why the first equality holds?


title: Finding the most relevant factors that lead to an event

text: I have been tasked with looking at a set of data with about 6000 records that each have 60 or so qualities associated with them (Let these be X1, X2, ... ) and determining what are the top 8 factors that determine whether a record will have a certain designation (there are three, Let's say, A, B, and C). Most of these X's only have two possible values so they are easy to deal with. However there is an X that has 6 values and an X that has 8 values which I have determined to be fairly important. For now, I am just interested in finding the factors amongst the data that lead a record to be event A. My approach so far has been to calculate all of the probabilities of each individual factor (A, B, C and all X's), calculate the probabilities of all of the X's given factor A, and to then use Bayes' theorem to calculate the probability of event A given the X. This gives me a feel for which factors seem to be contributing to event A but I am not sure how I should continue my analysis. I also delved into looking at intersections of events but I feel like that is a big time sink that yields little to no insight because that requires counting all of the intersections and when you want to look at 5+ factors, the amount of counting you have to do is ridiculous (For example, if you want to look at 5 factors, you need to find 2^5 different events because of the possible ways the events can line up. This causes problems especially when I am looking at the factors that have 6 or 8 different possibilities). I wrote some simple java classes to read in a tab delimited text file of the records and to do all of the counting and calculating for me for individual events and some intersections and I also have an excel spreadsheet that I set up first but that seemed to be fairly slow for performing the calculations. Besides adding specific code (or specific formulas in excel) to pick up each individual event, I am not sure how I can efficiently calculate these probabilities. And I'm not sure if calculating these conditional probabilities of the intersections will help me answer the overall question of what are the top 8 factors that contribute to event A. Any thoughts would be greatly appreciated.


title: Ordinal trends and finding concordant and discordant pairs

text: > A pair is concordant if the subject ranked higher on X also ranks higher on > Y. The pair is discordant if the subject ranking higher on X ranks lower on > Y. The pair is tied if the subjects have the same classification on X and/or > Y. Say I have a table (e.g., $n_{11}$ is the first cell: $i=1$, $j=1$) \begin{array}{rrrr} 1 &3 &10 &6 \\\ 2 &3 &10 &7 \\\ 1 &6 &14 &12 \\\ 0 &1 &9 &11 \end{array} Can you help me to understand how to select the concordant and discordant pairs? Here's my understanding for concordant:   First we select $n_{11}$, then $n_{22}$ is $3$ which is greater than $1$, (higher $Y$ and higher $X$). So we have $1(3+10+7+6+14+12+1+9+11)$, is that right? Then how about $n_{13}$? For $n_{13}$, we have $7$ ($n_{24}$) which is less than $n_{13}$ (lower $Y$). So we don't count this pair right? I am confused about discordant, how could you select the discordant pairs?


title: relation between stochastic independence and functional independence

text: We know if two variables are functionally independent then there will be stochastic independence as well. Can any one give me an example when the converse is not true i.e. X and Y are independent (stochastic) but not functionally independent. Thanks a lot.


title: Why not always use bootstrap CIs?

text: I was wondering how bootstrap CIs (and BCa in barticular) perform on normally- distributed data. There seems to be lots of work examining their performance on various types of distributions, but could not find anything on normally- distributed data. Since it seems an obvious thing to study first, I suppose the papers are just too old. I did some Monte Carlo simulations using the R boot package and found bootstrap CIs to be in agreement with exact CIs, although for small samples (N<20) they tend to be a bit liberal (smaller CIs). For large enough samples, they are essentially the same. This makes me wonder whether there is any good reason not to _always_ use bootstrapping. Given the difficulty of assessing whether a distribution is normal, and the many pitfalls behind this, it seems reasonable not to decide and report bootstrap CIs irrespective of the distribution. I understand the motivation for not using non-parametric tests systematically, since they have less power, but my simulations tell me this is not the case for bootstrap CIs. They are even smaller. A similar question that bugs me is why not always use the median as the measure of central tendency. People often recommend to use it to characterize non normally-distributed data, but since the median is the same as the mean for normally-distributed data, why make a distinction? It would seem quite beneficial if we could get rid of procedures for deciding whether a distribution is normal or not. I am very curious about your thoughts on these issues, and whether they have been discussed before. References would be highly appreciated. Thanks! Pierre


title: Lack of fit and Pure error

text: I don't understand the concepts of lack of fit error and pure error. What I know is: $\bullet$ _Lack of fit error:_ Error that occurs when the analysis omits one or more important terms or factors from the process model. $\bullet$ _Pure error:_ I occurs for repeated values of dependent variable, Y for a fixed value of independent variable, X. Can you please explain these two terms to me? If there is no repeated observations, does the pure error occur? If not, can't it possible to occur lack of fit error? Do they need to occur simultaneously ? Can't lack of fit error solely contribute to residual ? that is , $$\text{Residual Error=Lack of fit error + Pure Error}$$ If Pure Error=0, then can't be it $$\text{Residual Error=Lack of fit error}$$ Since we do the test $$F=\frac{\text{Mean Square due to Lack of fit}}{\text{Mean Square due to Pure error}}$$ It is good to have larger value of Pure error. Doesn't it imply our observation's are so heterogeneous. Isn't it good to have homogeneous units?


title: Meta-analysis of logistic regression

text: I would like to perform meta-analysis on studies that have a dichotomous outcome measure. I have the raw data for all studies I am including in my meta-analysis, so I can choose the best statistics to use. However, is not clear to me what statistics would be most appropriate for meta-analysis and I can't find any papers that discuss this issue. I've been using binary logistic regression in SPSS (but can switch to using R if necessary). In the SPSS model summary I get two measures of pseudo R-square -- Cox and Snell R square and Nagelkerke R square. These might be the best statistics to use in a meta-analysis. However, I've read (in this forum and elsewhere) that pseudo R square values aren't necessarily the most useful outputs generated by logistic regression. Perhaps it would be better to use the confidence intervals of the beta values for the predictors I am interested in? Or perhaps some other statistics? Any help would be much appreciated.


title: Understanding matlab code for converting raw input images into matlab format

text: f = fopen('train-images-idx3-ubyte','r'); [a,count] = fread(f,4,'int32'); g = fopen('train-labels-idx1-ubyte','r'); [l,count] = fread(g,2,'int32'); fprintf(1,'Starting to convert Training MNIST images (prints 60 dots)\n'); n = 1000; Df = cell(1,10); for d=0:9, Df{d+1} = fopen(['digit' num2str(d) '.ascii'],'w'); end; for i=1:60, fprintf('.'); rawimages = fread(f,28*28*n,'uchar'); rawlabels = fread(g,n,'uchar'); rawimages = reshape(rawimages,28*28,n); for j=1:n,               fprintf(Df{rawlabels(j)+1},'%3d ',rawimages(:,j));          fprintf(Df{rawlabels(j)+1},'\n');      end; end; fprintf(1,'\n'); for d=0:9, fclose(Df{d+1}); D = load(['digit' num2str(d) '.ascii'],'-ascii'); fprintf('%5d Digits of class %d\n',size(D,1),d); save(['digit' num2str(d) '.mat'],'D','-mat'); end; can somebody explain this code?


title: Showing characteristic function is absolutely continuous

text: A distribution has the characteristic function $\phi(t)$ = $(1-t^2/2)exp(-t^2/4), -\infty < t < \infty$ Show that the distribution is $absolutely$ continuous and write the density function of the distribution. Attempt: $\int_{-\infty}^{\infty}|(1-t^2/2)exp(-t^2/4)|dt $ = $(-2/t)(1-t^2/2)exp(-t^2/4)-2exp(-t^2/4)|_{-\infty}^{0}$ Similar result for $[0,\infty]$ since t is squared. I'm not quite sure I did the integration right, but if I can show that the absolute value of $\phi(t)$ is less than $\infty$, then the function is absolutely continuous.


title: Clustering data that has mixture of continuous and categorical variabes

text: I have data that represent some aspect of human behavior. I want to cluster it (unsupervised) into behavioral profiles of some sort. now, some of my variables are categorical (with 2 or more categories), and some are continuous (most are percentages). A few variables are even more complex in that one category has further continuous and the other one has no such additional data. My question is about how to go about categorize this data. What are the (common?) approaches dealing with it? I don't need code or anything, but rather some references or directions that will help me further understand how to deal with this challenge. If you know of `R` functions that facilitate such analysis, that would be great, but it's not necessary. thanks.


title: How do I analytically transform Wishart to inverse-Wishart?

text: How do I transform the probability density function of the Wishart distribution into that of the inverse-Wishart distribution? I'd like to know the whole process.


title: What is the test used in R in sm.regression to compare reference model with kernel regression?

text: I am doing bivariate kernel regression using the sm.regression function:   http://cran.r-project.org/web/packages/sm/sm.pdf There is an option to compare the nonparametric estimate with linear model. When requested, it returns p-value that the nonparametric estimate is not different from the linear (null hypothesis). Does anyone know what test is used to calculate the p-value? I can't find it anywhere. The answer might be in: Bowman, A.W. (2006). Comparing nonparametric surfaces. Statistical Modelling, 6, 279-299. Thank you!


title: Entropy of Cauchy (Lorentz) Distribution

text: Entropy is defined as $H$ = -$\int p(x)$ $log$ $p(x)$ $dx$ The Cauchy Distribution is defined as $f(x)$ = $\frac{\gamma}{\pi}$ $\frac{1}{\gamma^2 + x^2} $ I kindly ask to show the steps to calculate the Entropy of a Cauchy distribution, which is $log$($\gamma$) + $log$($4$$\tau$)


title: Why does boot.ci BCa fail on discrete data?

text: I've been using boot.ci to compute BCa interval estimates of medians but ran into problems when the data has too many identical values -- which can easily happen for ordinal data, but also discrete or low-precision interval data. If I add very small random numerical perturbations to the data in order to make all values different, BCa again works, and gives me intervals in agreement with the percentile method. I have two questions: 1) Is this workaround legitimate? 2) Can't the BCa code be modified to deal with these boundary cases? Although the percentile method works for these datasets, I would prefer not to have to switch methods depending on the error messages I get, and it's also much better to report a single method in a paper. BCa seems to be the best method to compute bootstrap confidence intervals, and is believed to work in most cases: Kirby KN, Gerlanc D. BootES: an R package for bootstrap confidence intervals on effect sizes. Behav Res Methods, 2013. Here is R code that illustrates this problem:               # Load package boot     library(boot)          # Make the code deterministic     set.seed(0)          # Define the median as the statistic of interest     samplemedian <- function(x, d) {return(median(x[d]))}          # Use many replicates     replicates <- 10000          # A problematic dataset     data1 <- c(1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5)          # The point estimate is 2     samplemedian(data1)          # boot.ci issues the following error:     #    Error in bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t0.o,  :      #      estimated adjustment 'w' is infinite     b <- boot(data1, samplemedian, R = replicates)     boot.ci(b, type = "bca")          # Add a negligible perturbation to the data     epsilon <- 1e-10     data1 <- data1 + runif(length(data1), -epsilon, epsilon)          # boot.ci now returns (2, 3) as the BCa confidence interval     b <- boot(data1, samplemedian, R = replicates)     boot.ci(b, type = "bca")          # Another problematic dataset     data2 <- c(1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5)          # The point estimate is 2     samplemedian(data2)          # boot.ci returns (1, 1) and issues the following warnings:     #    Warning : BCa Intervals used Extreme Quantiles     #    Some BCa intervals may be unstable     #    Warning message:     #    In norm.inter(t, adj.alpha) : extreme order statistics used as endpoints     b <- boot(data2, samplemedian, R = replicates)     boot.ci(b, type = "bca")          # Add a negligible perturbation to the data     epsilon <- 1e-10     data2 <- data2 + runif(length(data2), -epsilon, epsilon)          # boot.ci now returns (2, 3) as the BCa confidence interval     b <- boot(data2, samplemedian, R = replicates)     boot.ci(b, type = "bca")      Thanks in advance for your help! Pierre


title: Clusterings that can be caused by K-means

text: I have gotten the following question as a test question for my exam and I simply cannot understand the answer. A scatter plot of the data projected onto the first two principal components is shown below. We wish to examine if there exists some group structure in the data set. To do this, we have run the k-means algorithm with k = 2 using the Euclidean distance measure. The result of the k-means algorithm can vary between runs depending on the random initial conditions. We ran the algorithm several times and got some different clustering results. Only three of the four clusterings shown can be obtained by running the k-means algorithm on the data. Which one can not be obtained by k-means? (there's nothing special about the data) ![4 possible clusterings of data](http://i.stack.imgur.com/R9bIc.png) The correct answer is D. Can any of you explain why?


title: What are the implications of sampling to attain a given distribution in one variable?

text: Some while ago I attended a presentation of a case-control study for which enrolment was still in progress. The lecturer, a PhD student with basic statistical education, mentioned that while in the treatment group, age was already normally distributed, the corresponding distribution in the control group was still far from normal, and therefore they were actively searching for people with the "right" age to get that nice Gaussian curve. I have a strong feeling that proceeding in this fashion is at best problematic: Sampling this way is neither random nor independent, and age is normally distributed neither in the general population (from which the control group presumably is recruited) nor in patients with that disease (whose incidence is, of course, age-dependent). On the other hand, maybe it is not problematic at all, as it could be seen as simply matching cases to controls. I'm unsure of the implications for different analyses, in particular: Would one expect that inference on a treatment effect would be affected (age is a strong predictor of course of disease)?


title: Joint Cumulative Distribution Functions (with Marginal Properties)

text: I have encountered a joint cumulative distribution function with marginal properties. Although I have the answer, mine doesn't really match it and I don't seem to understand why. IMG: http://i.stack.imgur.com/0fwuN.jpg Following is my understanding, which seems to be in conflict of the resolutely answer. IMG: http://i.stack.imgur.com/QG7Fn.jpg Appreciate any advice please. PS: can't post image as I don't have enough rep points yet. :(( still learning the proper equations and syntax.


title: What technical language to describe the degree to which probabilities are likely to be modified by future data?

text: I'm trying to reason about something I call "estimate stability," and I'm hoping you can tell me whether there’s some relevant technical language, so that I can learn about it and then write a tutorial on the subject for Less Wrong. What do I mean by "estimate stability?" Consider these three different propositions:   1. We’re 50% sure that a coin (known to be fair) will land on heads.   2. We’re 50% sure that Matt will show up at the party.   3. We’re 50% sure that Strong AI will be invented by 2080. These estimates feel different. One reason they feel different is that the estimates have different degrees of "stability." In case (1) we don't expect to gain information that will change our probability estimate. But for cases (2) and (3), we may well come upon some information that causes us to adjust the estimate either up or down. So estimate (1) is more "stable," but I'm not sure how this should be quantified. Should I think of it in terms of running a Monte Carlo simulation of what future evidence might be, and looking at something like the variance of the distribution of the resulting estimates? What happens when it’s a whole probability distribution for e.g. the time Strong AI is invented? (Do you do calculate the stability of the probability density for every year, then average the result?) Here are some other considerations that would be useful to relate more formally to considerations of estimate stability:   * If we’re estimating some variable, having a narrow probability distribution (prior to future evidence with respect to which we’re trying to assess the stability) corresponds to having a lot of data. New data, in that case, would make less of a contribution in terms of changing the mean and reducing the variance.   * There are differences in model uncertainty between the three cases. I know what model to use when predicting a coin flip. My method of predicting whether Matt will show up at a party is shakier, but I have some idea of what I’m doing. With the Strong AI case, I don’t really have any good idea of what I’m doing. Presumably model uncertainty is related to estimate stability, because the more model uncertainty we have, the more we can change our estimate by reducing our model uncertainty.   * Another difference between the three cases is the degree to which our actions allow us to improve our estimates, increasing their stability. For example, we can reduce the uncertainty and increase the stability of our estimate about Matt by calling him, but we don’t really have any good ways to get better estimates of Strong AI timelines (other than by waiting).   * Value-of-information affects how we should deal with delay. Estimates that are unstable in the face of evidence we expect to get in the future seem to imply higher VoI. This creates a reason to accept delays in our actions. Or if we can easily gather information that will make our estimates more accurate and stable, that means we have more reason to pay the cost of gathering that information. If we expect to forget information, or expect our future selves not to take information into account, dynamic inconsistency becomes important. This is another reason why estimates might be unstable. One possible strategy here is to precommit to have our estimates regress to the mean. Thanks for any thoughts!


title: Multi-armed bandit algorithms in Java?

text: Multi-armed_bandit problem defenition from Wikipeda: "In probability theory, the multi-armed bandit problem (sometimes called the K-[1] or N-armed bandit problem) is the problem a gambler faces at a row of slot machines, sometimes known as "one-armed bandits", when deciding which machines to play, how many times to play each machine and in which order to play them. When played, each machine provides a random reward from a distribution specific to that machine. The objective of the gambler is to maximize the sum of rewards earned through a sequence of lever pulls" (http://en.wikipedia.org/wiki/Multi-armed_bandit) In other words Multi-armed bandit (MAB) problems are a class of sequential resource allocation problems concerned with allocating one or more resources among several alternative (competing) projects / 'arms'. Such problems have a conﬂict between making decisions (allocating resources) that yield high current rewards, versus making decisions that sacriﬁce current gains with the prospect of better future rewards. Any Java implementations for these algorithms?


title: A fixed effects regression in Stata that changes when a time invariant variable is included

text: I am running a fixed effects regression in Stata. If I understood it correctly, a variable that is fixed over time will be omitted from the regression. Hence, it should not influence the regression results. This is not the case, however. I ran a fixed effects regression with a time invariant variable. Stata omitted this variable from the regression. When I manually removed it and ran the regression again, the predicted outcome from another variable kept significant at the 5%, but changed sign! Any plausible explanation for this?


title: How do you interpret results from unit root tests?

text: I have to do some unit root tests for a project, I'm just unsure on how to interpret the data (which is what I have been asked to do). Here is one of my results:               dfuller Demand          Dickey-Fuller test for unit root                   Number of obs   =        50                            ---------- Interpolated Dickey-Fuller ---------                     Test         1% Critical       5% Critical     10% Critical                  Statistic         Value            Value           Value                  -------------------------------------------------------------------     Z(t)           -1.987         -3.580            -2.930          -2.600            -------------------------------------------------------------------               MacKinnon approximate p-value for Z(t) = 0.2924      What do I say about the critical values and the p-value results?


title: Autoregression model with a time trend term. Statistically valid?

text: Assume we have a time-series with a deterministic trend. I'm wondering if the following model is well specified as an AR model: $y_{t} = b_{0} + b_{1}t + b_{2}y_{t-1} + \epsilon_{t} \ \ \ \ \ \ $ Specifically, I'm interested in knowing if it's OK to put a time term and auto-regressive terms in the same model. If it is not, what would be the better way of specifying the model?


title: Which statistical method should I use to prove significant difference for my picture?

text: I would like to confirm my biological hypothesis with statistical approach. Red line on picture represents average frequency of some biological phenomenom for whole chromosome. Chromosome consists of milions of bases, which can be understood as positions. So that we have milions of positions total. Then I, using certain rule, pick up 100 positions. For these 100 positions I determine frequency of the same biological phenomenon (so that I have 100 isolated points, for every x there is y and they are connected by lines, as seen on picture). Because picture shows just these positions, red line does not represent average of picture (it represents average of whole chromosome). I want to show that freguency at my certain positions is significantly higher than average (so that the biological process I observe is really there). My data are from database,so my output is what I calculated using this database-there are no more measurements or data of this kind. y-axis represent counts, not ratios. Adjacent points on graph are neighbours on the chromosome, so one could expect that two neighbours on the chromosomes do experience same/similiar biological process. I am a computer scientist (with one course of statistics so I have an idea) - detailed howto would be appreciated. However, advice about which method to use may be fair enough :-) ![enter image description here](http://i.stack.imgur.com/8vtfF.jpg)


title: What if interaction wipes out my direct effects in regression?

text: In a regression, the interaction term wipes out both related direct effects. Do I drop the interaction or report the outcome? The interaction was not part of the original hypothesis.


title: How to interpret an insignificant main effect that became significant at the presence of interaction effects in moderated regression analysis?

text: > **Possible Duplicate:**   >  What if interaction wipes out my direct effects in regression? I did hierarchical regression analysis on my data due to having interaction effects in my research model. All the variables are continuous. ![enter image description here](http://i.stack.imgur.com/5BQB5.jpg) R2 increased from .695 in model1 (main effect only) to .734 in model2 (main &interaction effects)(sig. F change = .000). All the assumptions for the regression analysis have been met. I have two problems with the "coefficients" table, though: ![enter image description here](http://i.stack.imgur.com/AtYtX.jpg)   1. As u can see in the table, the insignificant beta value of **ZSC** in model1 became significant in model2! Is it ok? I'm confused! Which value should i consider to reject/accept the related hypothesis? B of model1 (which rejects the hypothesis) or 2 (which confirms it!!)?   2. Although the Beta value for **ZSC_X_CS** is significant, its positive sign is against the hypothesis! it's supposed to have a negative sign according to the literature & also logic! How should i treat this hypothesis? Accept? Reject? Partially accept?!!!


title: Interaction Suppresses Main Effects in Repeated Measures GLM?

text: > **Possible Duplicate:**   >  What if interaction wipes out my direct effects in regression? I have three hypotheses within my experiment, and interpreting the interaction in relation to the rest of the results is what I am having difficulty with.   1. The first is that the IV in the model (with two levels) will be related to one DV but not the other (time 1, but not time 2). A repeated measures GLM indicated that the IV was related to the _both_ DVs.   2. The second hypothesis was that a covariate will be related to one DV and not the other. Using a repeated measures GLM this hypothesis was supported.   3. Predicted an interaction between the IV, DV and a covariate. (A repeated measures GLM with IV (fixed factor), DVs (within factor), covariate). What I found was no effect of the covariate, no effect of the interaction, however my main effect for IV and DVs was gone (p=0.075) (hypothesis one). Now I know p=0.05 isn't everything, so I have two questions. How would I go about interpreting the fact that my main effect for hypothesis one disappeared IF:   1. you _don't_ consider p=0.05 a black and white cut off point?    2. you _do_ consider p=0.05 a black and white cut off point?  Does it change how I interpret the results from hypothesis one, or their meaning? Does it make it redundant?


title: Interaction term change the significance of the main effect variable

text: > **Possible Duplicate:**   >  What if interaction wipes out my direct effects in regression? I am just wondering why would adding an interaction term to the regression model change the main effect from significant to insignificant? I am running a poisson regression using a database (>40,000 subjects) I have 3 main explanatory variables, time (continous), gender and age group (4 groups). All the 3 expl. variables were significant before adding the interaction between time and age groups. after adding the inetraction, the time became insignificant (changed from p<0.0001 to 0.6) the age groups remained significant. The interaction term is significant only between time and one of the age groups. By the way, I tested the time centered and squared and it was not significant. Any suggestions?


title: What's the difference between mathematical statistics and statistics?

text: What's the difference of mathematical statistics and statistics? I've read this: > Statistics is the study of the collection, organization, analysis, and > interpretation of data. It deals with all aspects of this, including the > planning of data collection in terms of the design of surveys and > experiments. And this: > Mathematical statistics is the study of statistics from a mathematical > standpoint, using probability theory as well as other branches of > mathematics such as linear algebra and analysis. So what would be the difference beetween them? I can understand that the processes of collection may not be mathematical, but I guess that organization, analysis and interpretation are, am I missing something?